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Abstract 

Demonstration  of  a  strategy  to  perform  two-dimensional  diode  laser  tomography 
using  a  priori  knowledge  from  symmetry  arguments  and  computational  fluid  dynamic 
(CFD)  calculations  is  presented  for  a  flat  flame  burner.  The  strategy  uses  an  optimization 
technique  to  determine  flame  diameter  and  location  using  a  vector  quantization  approach. 
Next,  the  variance  in  a  training  set,  produced  from  CFD  calculations,  is  captured  using 
principal  components  analysis.  The  information  in  the  training  set  allows  interpolation 
between  beam  paths  resulting  in  temperature  and  density  maps.  Finally,  the  TDLAS 
temperature  and  density  maps  are  shown  to  agree  with  traditional  thermocouple 
measurements  of  the  flat  flame  burner  to  within  five  percent.  Preliminary  tomography 
results  using  the  same  strategy  are  then  presented  for  a  super-sonic  flow  through  an 
isolator,  a  length  of  wind  tunnel  necessary  to  keep  the  shock  wave  produced  by  the 
combustor  from  reaching  the  inlet. 
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I  Introduction 


Background 

Since  the  mid  1990s,  groups  have  used  tunable  diode  laser  absorption 
speetroseopy  (TDLAS)  to  determine  density,  pressure,  temperature,  and  velocity  of 
gases.[l]  TDLAS  uses  a  spectrally  narrow,  tunable  laser  source  to  scan  an  absorption 
feature  of  a  particular  gas.  The  ability  to  measure  gas  density,  pressure,  temperature,  and 
velocity  in  the  flow  of  a  supersonie  eombustion  ramjet  (scramjet)  engine  would  allow 
engine  performanee  to  be  monitored.  Optimizing  engine  performanee  in  a  scramjet 
presents  many  teehnieal  ehallenges  including  determination  of  mass  eapture  as  a  function 
of  vehicle  attitude  and  condition  in  the  inlet  of  the  engine;  ignition,  flame  holding,  and 
fuel  injeetion  in  the  eombustor;  and  maximizing  combustion  efficiency  and  thrust 
generation  in  the  exhaust  plane. [2]  Although  other  teehniques  exist  for  obtaining 
information  about  the  flow,  TDLAS  can  be  used  to  monitor  all  of  these  variables  in 
ground  tests  and  also  has  the  potential  to  be  miniaturized  for  flight  vehieles. 

Alternative  techniques  involving  probe  and  wall  measurements,  while  proven,  do 
not  necessarily  provide  a  eomplete  pieture  of  the  flow.[3]  Wall  measurements  are  limited 
to  pressure  and  temperature  readings  at  the  boundary  of  the  flow.  Computational  fluid 
dynamie  (CFD)  calculations  and  probe  measurements  show  that  flow  eonditions  vary 
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greatly  from  the  core  of  the  flow  to  the  extremities.  TDLAS  improves  upon  wall 
measurements  by  providing  an  average  measurement  across  the  entire  line  of  sight  thus 
providing  information  in  the  core  of  the  flow.  Probe  measurements  are  beneficial  in  that 
they  can  be  moved  or  built  in  an  array  to  investigate  conditions  throughout  the  flow. 
However,  their  presence  perturbs  the  flow  necessitating  a  correction  to  the  measurement. 
TLDAS  does  not  have  this  drawback  because  passing  a  laser  through  the  flow  does  not 
disrupt  the  flow. 

TDLAS,  therefore,  provides  another  avenue  for  collecting  information  about  the 
flow.  By  passing  a  beam  through  the  flow,  a  line  of  sight  absorption  measurement  is 
possible.  By  carefully  choosing  the  laser  wavelengths  to  correspond  to  species- specific 
absorption  features,  of  water  in  the  case  of  this  study,  a  variety  of  flow  properties  can  be 
determined.  Three  lasers  with  nominal  center  wavelengths  of  1391,  1392,  and  1397  nm 
were  used  along  each  line  of  sight  in  this  experiment.  The  shape  and  intensity  of  the 
absorption  features  in  each  of  the  three  spectra  yielded  information  regarding  temperature 
and  pressure  as  well  as  the  density  and  velocity  of  the  species  in  the  probed  area  of  the 
flow.  [3] 

Line  of  sight  measurements  are  limited,  however,  in  that  the  information  gained  is 
averaged  along  the  path  of  the  beam.  With  a  single  beam  path,  therefore,  no  spatial 
information  about  the  flow  can  be  obtained.  One  solution  to  gathering  spatially  resolved 
information  of  the  flow  from  the  line  of  sight  measurements  is  to  add  more  beam  paths. 
Tomography,  similar  to  that  used  in  an  MRI  or  CAT  scan,  is  then  performed  by 
combining  the  information  of  multiple  beam  paths.  Unfortunately,  methods  like  these 
generally  require  hundreds  of  beam  paths  to  gain  the  resolution  necessary  to  extract 
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contour  information  -  a  solution  which  is  not  feasible  in  flight  environments  or  even  most 
ground  tests.  To  overcome  this  problem,  reduction  in  the  number  of  beam  paths 
necessary  to  perform  tomography  can  be  accomplished  through  the  use  of  contour 
information  gained  from  CFD  computations  or  planar  laser-induced  fluorescence  (PLIF) 
experiments.  By  assuming  the  model  gained  from  CFD  or  PLIF  provides  a  good 
approximation  of  the  relative  variation  in  the  flow  properties,  the  model  can  be  fit  to  the 
TDLAS  line  of  sight  information  through  a  least  squares  fitting  method,  allowing 
intelligent  interpolation  between  the  beam  paths.  Through  this  process  of  two- 
dimensional  tomography,  spatially  and  temporally  resolved,  quantitative  information 
about  the  flow  can  be  obtained.  [3,  4] 

Problem  Statement 

The  overall  purpose  of  this  project  was  to  investigate  the  possibility  of  reducing 
the  number  of  lines  of  sight  required  to  perform  tomography  using  a  priori  knowledge 
gained  from  symmetry  arguments  or  CFD  calculations.  To  fulfill  this  purpose,  the 
temperature  and  water  concentration  of  a  flame  provided  by  a  flat  flame,  McKenna  type 
burner  were  measured  six  millimeters  above  the  surface  of  the  burner  through  the 
traditional  method  of  thermocouples  as  well  as  a  pressure  and  relative  humidity  probe. 
The  flat  flame  burner  consists  of  a  circular,  porous  copper  region,  six  centimeters  in 
diameter.  The  gas  mixture,  methane  and  air,  is  distributed  evenly  over  this  region  and, 
when  lit,  provides  a  uniform,  or  flat,  temperature  profile  across  the  area  of  the  porous 
region.  The  results  from  the  thermocouple  measurements  were  compared  to  the  tunable 
diode  laser  absorption  tomography  results  performed  using  an  analytic  two-zone  model, 
an  optimized  computational  two-zone  model,  and  finally  an  optimized  computational 
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CFD  model  of  the  flat  flame  burner.  The  two-zone  model,  discussed  in  detail  in  Chapter 
IV,  assumes  two  temperature  regions  for  the  flat  flame  burner  experiment:  an  elevated 
temperature  within  the  flame  of  the  burner  surrounded  by  room  temperature  air.  The 
agreement  on  temperature  and  species  concentration  of  the  two  methods,  tunable  diode 
laser  absorption  tomography  and  thermocouple  measurements,  provides  a  proof  of 
concept  for  the  tomography  strategy  which  was  employed.  This  same  strategy  was  then 
applied  to  an  experiment  to  perform  tunable  diode  laser  absorption  tomography  (TDLAT) 
on  a  vitiated  supersonic  flow  through  an  isolator.  Preliminary  results  are  discussed  in 
Chapter  V. 

The  two-zone  model  was  used  in  the  implementation  of  two  methods  to  determine 
temperature  and  species  concentration.  The  first  method,  called  the  two-zone  analytic 
method,  uses  TDLAS  data  from  seven  beam  paths,  the  positions  of  the  lines  of  sight 
relative  to  the  burner,  and  two  assumptions  about  the  burner:  the  flat  flame  burner 
provides  a  flat  temperature  profile,  and  the  flame  is  circular.  Employing  some  algebra 
and  geometry,  detailed  in  Chapter  IV,  the  temperature  in  the  center  of  the  flame  was 
found  to  be  within  10  K  of  the  thermocouple  measurements.  The  second  method,  called 
the  two-zone  computational  method,  used  the  same  data,  beam  path  measurements,  and 
assumptions  about  the  burner  to  obtain  the  same  temperature  and  water  concentration 
results  computationally.  In  the  computational  method,  the  two-zone  model  provided  a 
basis  image  which  was  scaled  to  fit  the  TDLAS  data  using  a  least  squares  method  also 
detailed  in  Chapter  IV.  These  two  methods  are  intended  to  step  the  reader  through  a 
basic  application  of  tomography  using  a  basis  derived  from  a  priori  knowledge.  Both  of 
these  methods  agreed  with  the  thermocouple  temperature  measurements  and  species 


4 


concentration  calculations  to  within  five  percent  for  the  interior  of  the  flame  and  the 
surrounding  room  air,  but  because  they  use  the  two-zone  model,  they  do  not  provide  any 
information  at  the  boundary  layer  of  the  flame.  The  two-zone  model  assumes  an 
instantaneous  temperature  ehange  at  the  boundary  between  the  surrounding  room 
temperature  gas  and  the  hot  eircular  flame.  To  obtain  better  agreement  with  the 
thermocouple  measurements  in  this  region,  either  more  beam  paths  along  with  a  model 
with  more  temperature  zones  or  a  better  understanding  of  how  to  interpolate  between  the 
beam  paths  was  necessary.  CFD  ealeulations  were  used  to  provide  this  information. 

The  CFD  model  provided  the  information  to  allow  interpolation  between  lines  of 
sight,  thus  obtaining  spatially  resolved  temperature  and  speeies  eoneentration  maps  of  the 
flat  flame  burner.  Although  the  temperature  and  diameter  of  the  CFD  ealeulations  did 
not  agree  with  thermoeouple  measurements,  the  shape  of  the  temperature  profile  in  the 
boundary  region  showed  good  agreement  with  the  thermocouple  data.  By  sealing  the 
CFD  model  to  represent  a  range  of  peak  temperatures  and  peak  speeies  eoneentrations  of 
the  burner,  a  training  set  was  developed  for  each  of  the  water  absorption  features  used  for 
temperature  and  species  eoneentration  determination.  The  training  sets  for  each  of  the 
absorption  features  consisted  of  a  series  of  maps  of  the  frequeney  integrated  absorption 
coeffieients  for  a  range  of  sealed  temperatures  and  species  eoneentration. 

Two  algorithms  were  proposed  to  find  the  optimum  basis  set  from  the  training  set: 
vector  quantization  (VQ)  and  principal  components  analysis  (PCA).[5,  6]  Both  of  these 
algorithms  can  be  found  in  Appendix  C  in  the  programs  Create  VQbasis  and 
Create_basisv2,  respectively.  While  PCA  performed  well  in  the  ereation  of  basis  sets, 
VQ  showed  limited  utility  for  the  flat  flame  burner.  Instead,  VQ  was  used  to  optimize 
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the  fits  to  the  two  models,  to  aid  in  providing  averages  in  time-resolved  oscillations  in  the 
isolator  experiment,  and  to  create  a  basis  for  the  near  real-time  tomography  performed 
during  data  collection  for  the  isolator  experiment.  PCA  was  used  to  create  the  basis  sets, 
capturing  variance  in  each  training  set  and  resulting  in  a  basis  set  which  could  be  reduced 
in  size. [5]  A  least  squares  approach  was  used  to  fit  the  CFD  basis  sets  to  the  TDLAS 
measurements.  The  basis  sets  were  then  added  in  linear  combination  resulting  in  a 
tomographic  map  which  agreed  with  the  TDLAS  data  from  the  multiple  lines  of  sight. 

The  PCA  was  performed  because  the  CFD  models  covering  all  of  the  possible 
states  could  be  very  large,  thus  requiring  many  lines  of  sight  to  make  the  fit  well  defined. 
To  reduce  the  number  of  lines  of  sight  necessary,  PCA  was  used  on  the  training  sets  to 
create  much  smaller  basis  sets  which  still  represented  the  majority  of  possible  conditions 
of  the  flow  yet  maintained  a  well  posed  mathematical  problem  with  fewer  lines  of  sight. 
The  compression  of  the  training  sets  was  possible  due  to  the  high  degree  of  redundancy 
in  the  varied  conditions  of  the  CFD  models.  The  research  questions  under  investigation 
for  this  thesis  were  how  to  create  the  best  possible  basis  set  from  the  CFD  models  and 
how  to  obtain  the  best  possible  fit  of  the  measured  TDLAS  data  to  the  basis  set  while  still 
maintaining  a  well  conditioned  problem. 

Scope 

This  thesis  investigated  a  strategy  to  perform  tomography  of  a  flat  flame  burner 
and  vitiated  super-sonic  flow  using  TDLAS  measurements  from  multiple  lines  of  sight 
resulting  in  quantitative  temperature,  density,  and  pressure  maps.  Temperature  and 
species  concentration  maps  of  the  flat  flame  burner,  obtained  using  tomography  were 
compared  to  a  thermocouple  temperature  proflle  and  calculated  species  concentrations  of 
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the  flame.  This  comparison  provided  verification  for  the  fit  of  the  TDLAS  measurements 
to  the  two-zone  and  CFD  models.  The  favorable  comparisons  of  both  the  temperature 
and  species  concentration  obtained  using  tomography  to  those  obtained  using 
thermocouples  lend  added  credibility  to  the  optimized  fitting  routine  applied  to  the  super¬ 
sonic  flow  data  by  Lindstrom.[3]  Preliminary  results  for  the  super-sonic  flow  are 
presented  in  chapter  6. 

Objectives 

In  this  thesis,  the  following  objectives  were  accomplished: 

i.  Verify  the  utility  of  TDLAS  for  the  flat  flame  burner  by  comparing  with 
thermocouple  measurements. 

ii.  Investigate  the  number  of  beam  paths  required  to  achieve  accurate 
tomographic  temperature  and  density  maps  for  the  flat  flame  burner. 

iii.  Evaluate  the  performance  of  PC  A  in  reducing  the  dimensionality  of  training 
sets  and  enabling  tomography  with  a  limited  number  of  beam  paths. 

iv.  Provide  a  proof  of  concept  for  the  application  of  TDLAT  performed  with 
limited  optical  access  in  more  complicated  flows. 

Approach/Methodology 

A  temperature  profile  of  a  flat  flame  burner  measured  with  thermocouples  was 
used  as  a  comparison  for  the  TDLAT.  Tomography  was  performed  using  both  the  two- 
zone  model  as  well  as  the  CFD  model  using  the  data  from  7,  63,  and  1 19  lines  of  sight  for 
each  model.  This  provided  a  proof  of  concept  for  the  strategy  developed  to  perform  two- 


7 


dimensional  tomography  of  the  flat  flame  burner.  The  proof  of  concept  was  used  to 
provide  confidence  in  the  tomography  results  performed  on  the  isolator  data. 

Equipment  and  Materials 

The  flat  flame  burner  experiments  were  performed  in  the  Propulsion  Directorate 
of  the  Air  Force  Research  Laboratory  at  Wright-Patterson  Air  Force  Base.  The  first 
experiment  measured  the  temperature  profile  of  a  flat  flame,  McKenna  type  burner  using 
an  array  of  three  type  B  thermocouples  of  varying  bead  and  wire  diameters  and  two  type 
K  thermocouples.  The  second  burner  experiment  was  performed  with  an  array  of  tunable 
diode  lasers  and  the  required  data  acquisition  system  to  measure  water  absorption 
features  which  were  used  to  obtain  a  similar  temperature  profile. 

The  isolator  experiments  were  performed  in  Research  Cell  1 8  at  Wright-Patterson 
Air  Force  Base.  These  experiments  used  the  same  tunable  diode  lasers  and  data 
acquisition  system  as  the  TDLAS  burner  experiment  to  measure  water  absorption  in  the 
flow  through  the  isolator  portion  of  a  supersonic  wind  tunnel. 

Data  processing  and  tomography  were  performed  on  a  standard  desktop  computer 
with  a  2.39  GHz  processor  and  2.0  GB  of  RAM.  Algorithm  development  was 
accomplished  primarily  in  Matlab,  though  the  Igor  curve  fitting  suite  was  used  to  process 
the  raw  data  from  the  tunable  diode  lasers. [7] 

Preview 

The  remainder  of  this  thesis  is  organized  into  five  chapters.  Chapter  2  provides 
the  necessary  spectroscopic  background  and  theory  to  obtain  a  tomographic  temperature 
map  from  raw  TDLAS  measurements.  Chapters  3,  4,  and  5  provide  the  background, 
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experimental  setup,  and  results  for  each  of  the  three  experiments:  the  flat  flame 
thermocouple  experiment,  the  flat  flame  TDLAT  experiment,  and  the  isolator  TDLAT 
experiment.  Chapter  3  explains  the  traditional  thermocouple  methods  used  to  obtain  a 
temperature  profde  for  the  flat  flame  burner.  The  results  of  the  thermocouple 
measurements  provide  a  point  of  comparison  for  the  remaining  TDLAS  measurements 
and  tomography.  Chapter  4  includes  the  results  for  both  the  thermocouple  and  TDLAS 
experiments  on  the  flat  flame  burner  along  with  the  comparison  of  the  two.  This 
comparison  validates  the  strategy  implemented  to  perform  two-dimensional  tomography 
on  a  flow.  Chapter  5  provides  the  preliminary  results  of  the  isolator  experiment  which 
employed  the  same  strategy  to  perform  two-dimensional  tomography  as  the  flat  flame 
burner.  Using  only  line  of  sight  TDLAS  measurements,  the  CFD  was  shown  to  be 
missing  the  first  reflection  in  the  shock  train.  Due  to  this  inaccuracy  in  the  CFD 
calculations  for  the  isolator,  optimization  techniques  used  in  the  flat  flame  experiment  to 
determine  diameter  and  center  location  were  modified  to  adjust  the  shock  shape,  location, 
and  wall  thickness  in  the  isolator  CFD  model.  The  resulting  fits  to  the  isolator  TDLAS 
data  are  presented.  The  thesis  concludes  with  Chapter  6  where  a  summary,  overview  of 
future  studies,  and  recommendations  for  future  research  are  discussed. 
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II  Theory 


Introduction 

This  chapter  develops  the  theory  behind  TDLAS  and  the  tomography  performed 
using  TDLAS  measurements.  It  begins  with  a  description  of  the  water  moleeule 
vibrational  and  rotational  energy  levels  that  provide  the  transitions  whieh  were  used  for 
TDLAS.  Next,  Beer’s  law  is  presented.  It  governs  the  absorption  of  light  by  the  water 
transitions.  The  ehapter  then  provides  an  overview  of  the  speetroseopie  background 
material  necessary  to  understand  TDLAS.  Finally,  a  discussion  of  tomography  as  it  is 
applied  to  this  thesis  is  presented. 

The  Water  Molecule 

The  burner  and  isolator  experiments  rely  on  absorption  speetra  which  are  the 
result  of  the  absorption  of  laser  light  by  water  molecules  in  a  flow.  As  the  laser  is  tuned 
through  a  ffequeney  range,  transitions  in  the  water  molecule  attenuate  the  light  aceording 
to  the  temperature  and  concentration  of  the  water  molecules.  Water  was  seleeted  as  the 
speeies  of  measurement  due  to  its  abundance  of  transitions  in  the  range  of  wavelengths 
where  eommercial  telecommunications  lasers  have  been  developed.  From 
thermodynamics,  molecules  have  3N  degrees  of  freedom  where  N  is  the  number  of  atoms 
in  the  moleeule.  For  the  water  molecule,  three  of  its  nine  degrees  of  freedom  are 
translational  and  three  are  rotational,  leaving  three  vibrational  modes.  Eaeh  vibrational 
state  is  eomposed  of  many  rotational  states.  The  figure  below,  adapted  from  Bemath, 
shows  the  vibrational  modes  of  the  water  moleeule.  [8] 
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Figure  1.  The  normal  modes  of  vibration  of  H2O 

While  many  water  vibrational  transitions  exist  within  the  wavelength  range  of  the  lasers 

used,  the  transitions  used  in  this  research  occurred  between  rotational  levels  in  the  000 

and  101  vibrational  transition.  The  three  numbers  in  the  vibrational  level  labeling 

scheme  refer  to  the  quanta  of  energy  in  each  of  the  vibrational  modes,  symmetrie  stretch, 

bend,  and  antisymmetric  stretch  respectively.  [8]  The  transition  from  000  to  101, 

therefore,  indieates  that  the  absorption  features  used  were  the  result  of  ineident  light 

exciting  water  moleeules  from  their  vibrational  ground  state  to  an  exeited  state  with  one 

quanta  of  symmetrie  streteh  and  one  quanta  of  antisymmetric  stretch. 

Specific  rotational  transitions  of  the  000  to  101  vibrational  transition  were 
seleeted  based  on  the  population  density,  obtained  using  Boltzmann  statisties,  of  the 
ground  vibrational  state.  For  the  molecule  to  absorb  light  at  a  transition,  it  must  have  a 
population  in  the  initial  state.  Thus,  transitions  were  seleeted  based  on  the  population  of 
the  initial  state  over  the  expected  temperature  range  of  the  experiment. 

In  rotational  spectroscopy,  energy  levels  are  referred  to  by  their  total  angular 
momentum,  J,  and  the  projeetion  of  the  total  angular  momentum  on  the  molecular  axis  of 
highest  symmetry,  z.  As  shown  below  in  the  figure  adapted  from  Bemath,  water  is  a 
planar  moleeule  where  the  highest  degree  of  rotational  symmetry  is  about  the  180° 
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rotationally  symmetric  axis  through  the  oxygen  atom  in  the  plane  of  the  molecule. [8] 
For  a  planar  molecule,  the  x-axis  is  chosen  to  be  out  of  the  plane. 


z 

y 

!  f 

I 

I 

I 

Figure  2.  The  water  molecule 

The  origin,  which  defines  the  distances  f,  g,  and  h,  is  chosen  to  be  the  center  of  mass  of 
the  molecule.  The  moments  of  inertia  for  the  rigid  rotor  model  of  the  water  molecule  are 
therefore: 

ly  =moh^ +2mfjg^  .  (1) 

I^=I,+Iy=m^h^+2mAg"+f") 

Axes  are  also  labeled  according  to  the  magnitude  of  their  respective  moments  of  inertia 
such  that 

(2) 

In  the  case  of  H2O  where  r  =  .0958  A,  0  =  104.5°,  mn  =  1.00  atomic  mass  unit  (u),  mo  = 
16.00  u,  f  =  0.0757  A,  g  =  0.5213  A,  and  h  =  0.0652  A,  the  moments  of  inertia  are 

7^=1.148  uk\^I,) 

7^=0.6115  (3) 

7^  =1.760  uk\=Ic) 
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The  water  molecule  is  thus  an  asymmetric  top  because  none  of  the  moments  of  inertia  are 
equal.  Asymmetric  tops  provide  a  large  number  of  transitions  because  the  dipole 
moments  along  each  of  the  axes,  a,  b,  and  c,  are  non- vanishing.  [8] 

Once  the  axes  are  defined,  the  rotational  energy  levels  and  selection  rules  for 
linear  and  symmetric  top  molecules  can  be  developed  using  the  Hamiltonian  operator  for 
the  rigid  rotor  in  the  Schrodinger  equation.  In  the  case  of  asymmetric  tops,  however, 
there  is  no  general  analytical  solution  to  the  Schrodinger  equation.  Instead,  solutions  are 
usually  found  numerically. 

For  this  research,  HITRAN  was  used  to  provide  the  transition  parameters 
including  initial  state  energy,  transition  energy,  and  line  intensity.  HITRAN  is  a 
collection  of  spectroscopic  parameters  created  by  the  Air  Force  Cambridge  Research 
Laboratories. [9]  The  line  intensity  of  a  transition.  S',  is  shown  below  in  SI  units  of 
m  s  .  It  is  proportional  to  the  evaluation  of  the  dipole  moment  operator,  Sjo,  in  the 
Schrodinger  equation.  It  combines  the  probability  of  a  transition  with  thermodynamic 
equilibrium  Boltzmann  statistics  given  by:  [8] 


^10 

3€Qhcq(T) 


-Eq  / kT  /-I  -hv^n  I kT  \ 

e  “  (1-e  )5io. 


Here,  vw  is  the  transition  frequency  between  rotational  levels  of  the  ground  vibrational 
state,  0,  and  the  upper  vibrational  state,  1,  A  is  Planck’s  constant,  c  is  the  speed  of  light, 
q(T)  is  the  partition  function,  Eg  is  the  energy  of  the  initial  state,  or  term  energy,  k  is 
Boltzmann’s  constant,  T  is  the  temperature,  and  Sjo  is  the  line  strength  given  by  the 
dipole  moment  operator.  The  line  intensity,  S',  can  be  related  to  the  cross  section  of  the 
transition  as  well  as  the  absorption  coefficient  and  the  Einstein  coefficients  as  shown  in 
the  next  section. [8] 
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Beer’s  Law 


When  light  passes  through  a  gas,  absorption  or  stimulated  emission  may  occur 
when  the  energy  of  the  incoming  photons  matches  the  energy  gap  between  two  levels.  A 
transition  from  a  rotational  state,  of  the  initial  vibrational  state  to  a  rotational  state  of  the 
excited  vibrational  state  is  shown  in  the  figure  below,  adapted  from  Bemath.[8]  The 
energy  of  the  incoming  photon,  E,  is  equal  to  he  divided  by  the  wavelength,  X,  as  shown 
in  Figure  3. 

1 

0 

Figure  3.  A  transition  of  the  water  molecule  from  the  initial  state,  0  to  the  upper  state,  1 

Beer’s  law  describes  the  intensity,  /,  of  light  after  it  passes  through  a  distance,  /,  of  gas. 
Given  an  incoming  radiation  density,  p,  and  the  Einstein  coefficients  for  stimulated 
emission  and  absorption,  and  ,  the  rate  equation  for  the  upper  state  is: 

(5) 

at 

where  Ni  and  No  are  the  populations,  given  in  molecules-cm'  ,  of  the  excited  and  initial 
state  respectively.  The  incoming  radiation  density,  p  =  hv\oF/c,  has  units  of  kg-s'  -m' , 
where  h  is  Planck’s  constant  in  J-s,  vjo  is  the  transition  frequency  in  Hz,  F  is  the 
incoming  flux  in  photons  m'  -s'  ,  and  c  is  the  speed  of  light  in  m-s'  .  Note  that  the 
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Einstein  coefficients  for  absorption  and  stimulated  emission,  and  ,  both  having 
units  of  m-s-kg  '  from  the  following  equation: 

where  go  =  2Jo  +  1  and  gi  =  2 J/  +  1  are  the  degeneracies  of  state  0  and  1  with  total 
angular  momentum  Jo  and  Ji,  is  the  permittivity  of  free  space  in  C^-N'^-m'^,  h  is 

Planck’s  constant  in  J-s,  g(v-vio)  is  the  lineshape  function  with  units  of  seconds,  and  Sjo  is 
the  line  strength  in  C^-m^: 

^o=ZZI(-^.^J>'koA^o)r-  (7) 

Ml  Mo 

In  Equation  (7),  Jo  and  Ji  are  the  total  angular  momentum  of  a  molecule  in  state  0  or  1 
and  Mo  and  Mj  are  the  projection  of  Jo  and  Ji  onto  the  lab  reference  frame  for  those 
states.  Rearranging  Equation  (5)  and  using  Equation  (6)  yields: 


=  (8) 

dt  gi 

Replacing  the  Einstein  absorption  coefficient  with  its  definition  from  Equation  (6) 
produces: 


dN^ 

dt 


P  7^^  S„g(v-vJ(N,-^N0^ 


Adding  the  definition  ofp  =  hvioF/c  generates: 


dN^ 

dt 


'is^hcg^ 


^N,). 

Sx 


Then,  the  cross  section,  a,  with  units  of  m^  is  defined  as  follows: 


(7  = 


In  iSjQ  ^ 
3£^hc  go 


(9) 


(10) 


(11) 
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with  Sio  in  C  -m  ,  vio  in  Hz,  g(v-vio)  in  seconds,  So  in  C  -N'  -m'  ,  h  in  J  s,  and  c  in  m-s'  . 

Substituting  the  cross  section  from  Equation  (11)  into  Equation  (10)  results  in  the  form 
below:  [8,  10] 

^  =  oF(N,-^N,).  (12) 

at  gi 

Once  the  rate  equation  is  defined  as  shown  in  Equation  (12),  the  derivation  of 
Beer’s  law  follows  by  considering  the  change  in  incident  flux  caused  by  an  area  of  one 
square  meter  with  thickness  dx  as  shown  in  the  figure  below,  adapted  from  Bernath.[8] 


Figure  4.  A  volume  of  gas 


In  this  case,  the  rate  of  change  of  flux  over  the  distance  dx  is: 


dF  =  -ctF{N^  -^N,)dx.  (13) 

The  result  after  integration  is: 

(  f'\  (  t'\ 

In  —  =ln  -  =-c7{N, (14) 

voj  V-^oj  §i 

where  /  is  the  path  length  through  the  gas,  F  =  hv/I,  and  Fg  =  hv/Ig.  Exponentiation  of 
each  side  of  Equation  (14)  gives  the  common  form  of  Beer’s  law: 

I  =  Ige  .  (15) 
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Finally,  the  HITRAN  line  intensity  can  be  related  to  the  line  strength  in  Equation 
(7)  by  the  following  formula: 

3S(fhcq(T) 

Equation  (16)  gives  the  HITRAN  line  intensity  in  SI  units  of  m^-s'\  The  HITRAN 
literature,  however,  uses  units  of  cm‘V(molecules  cm'^).  This  form  is  shown  below:  [9, 


^  X 10-'' . 

3hcq(T)  ^10 


where  v  has  units  of  cm' 


^,0=—,  (18) 
go 

and  S]o  from  Equation  (7)  is  given  in  Debye^  instead  of  the  SI  units  of  C^-m^.  Assuming 
thermodynamic  equilibrium,  Maxwell-Boltzmann  statistics  are  used  to  eliminate  and 

Aj  from  Equation  (15)  with  the  following  relationships: 


-'^l  _  8 1  ^-hVi„/kT 

^0  8o 


In  Equation  (19),  Nj  is  the  total  species  concentration  in  molecules-cm''  and  q(7’)  is  the 
partition  function.  Using  Equation  (19)  to  rewrite  Equation  (15)  results  in  the  form  of 
Beer’s  law  shown  here: 


\IJ  3s^hcq{T) 


Using  the  definition  of  S'  from  Equation  (17)  yields: 
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=  S'g{v-v,,)N,l. 


(21) 


Integrating  over  a  single  absorption  line  eliminates  the  line  shape  function  resulting  in  the 
final  useful  form: 


/•  00  ^00 

[  A(v)dv  =  [  -  In 

J  —00  J  —00 


V^oy 


\dv  =  S'Nj,l  =  a-l 


(22) 


where  A(v)  (unitless)  is  the  natural  logarithm  of  the  ratio  of  the  intensities  of  the  incident 
light  to  exiting  light  measured  by  the  tunable  diode  laser  along  its  line  of  sight  and  a  is 
the  frequency  integrated  absorption  coefficient.  If  v  is  reported  in  wavenumbers,  then  the 
integral  of  A  has  units  of  wavenumbers,  S'  has  the  HITRAN  units  of  cm  V(molecuIes  cm' 
^),  N  is  given  in  molecules/cm^,  /  is  given  in  cm,  and  a  has  units  of  cm'^.  A  more 
detailed  development  stepping  from  Equation  (15)  to  Equation  (22)  is  given  in  Appendix 
A.  The  plot  of  A(v)  versus  v  is  the  absorption  spectrum  for  a  laser  along  a  line  of  sight.[8, 
12] 


Tunable  Diode  Laser  Absorption  Spectroscopy 

TDLAS  uses  a  tunable  diode  laser  to  measure  the  absorption  of  water  transitions 
over  a  path  length.  Beer’s  law  can  then  be  used  to  extract  temperature  and  density  from 
the  frequency  integrated  absorption  coefficients  of  two  or  more  transitions.  A  tunable 
diode  laser  consists  of  a  cavity  in  a  semiconductor  material  where  the  intensity  is 
controlled  by  the  input  current.  As  the  current  to  the  tunable  diode  laser  is  ramped,  the 
change  in  temperature  alters  the  cavity  length  of  the  laser.  Thus,  current  modulation 
permits  the  diode  to  be  tuned  across  a  wavelength  range.  The  laser  output  is  sent  through 
the  gas  of  interest  and  collected  on  the  other  side.  The  natural  logarithm  of  the  reference 
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beam  divided  by  the  spectrum  recorded  at  the  detector  gives  the  absorption  spectrum  for 
a  line  of  sight,  A(vJ.  Integrating  this  absorption  spectrum  over  frequency  provides  the 
integrated  absorption  along  the  path,  fjfvjdv,  as  seen  in  Equation  (22). 

Integrating  the  absorption  spectra  from  multiple  transitions  along  a  single  line  of 
sight  allows  an  average  temperature  to  be  determined.  This  was  accomplished  by 
scanning  three  lasers,  each  in  their  own  wavelength  range,  to  measure  multiple  absorption 
features.  As  discussed  at  the  beginning  of  the  chapter,  telecommunications  diode  lasers 
were  used  in  this  experiment,  so  the  probed  wavelength  regions  (approximately  1391  nm 
to  1397  nm)  fell  within  the  larger  wavelength  range  of  the  telecommunications  diode 
lasers.  Another  consideration  when  choosing  the  specific  nominal  wavelength  of  each 
diode  laser  was  the  term  energies  (the  energy  of  the  initial  state)  of  the  transitions  found 
around  the  nominal  wavelength.  In  thermal  equilibrium,  the  population,  and  therefore 
ability  to  absorb  incoming  photons,  of  the  initial  state  of  a  transition  depends  on 
temperature.  Thus,  if  the  term  energy  of  a  transition  is  too  high  to  be  filled  thermally  at 
the  temperatures  seen  in  an  experiment,  the  transition  will  not  provide  enough  absorption 
to  measure  with  TDLAS.  Likewise,  if  temperatures  in  an  experiment  are  high,  transitions 
with  low  term  energies  may  already  have  a  reduced  initial  state  population.  Even  higher 
temperatures  may  not  cause  a  significant  reduction  in  initial  state  population  and 
therefore  the  absorption  through  these  regions  would  not  be  changed.  Thus,  transitions 
were  chosen  which  had  initial  state  populations  that  produce  measurable  differences  in 
absorption  when  they  encounter  different  temperatures  in  the  temperature  ranges 
provided  by  the  experiments. 
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Equation  (24)  can  then  be  solved  for  the  temperature,  T,  using  the  HITRAN  line 
intensities  and  measured  frequency  integrated  absorptions. [12]  This  method  can  be 
performed  on  each  line  of  sight,  before  tomography,  using  an  analytic  two-ring  model 
approach,  or  after  tomography  for  each  individual  pixel.  Both  of  these  models  are 
discussed  further  in  the  TDLAS  Burner  results  section  in  Chapter  IV. 

The  second  method  of  temperature  determination  plots  two  or  more  integrated 
absorptions  against  the  term  energy  of  each  transition  in  what  is  called  a  Boltzmann  plot. 
Substituting  Equation  (17)  into  Equation  (22)  and  rearranging  results  in:  [3,  12] 
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The  Boltzmann  plot  is  generated  by  plaeeing  the  left  hand  side  of  Equation  (26)  on  the 
ordinate  and  the  term  energy,  Eg,  on  the  abscissa.  Omitting  the  stimulated  emission  term, 
ln[l  -  causes  at  most  only  a  few  percent  error  at  2500  K  and  was  therefore 

generally  omitted  from  the  calculation.  A  linear  fit  was  accomplished  for  the  integrated 
absorptions  where  the  slope  of  the  fit  is  equal  to  -  XlkT.  Again  this  method  was  performed 
analytically  for  each  laser  line  of  sight  and  pixel  by  pixel  for  the  absorption  tomographic 
maps.  Using  the  two-zone  analytic  model,  the  temperature  from  the  Boltzmann  plot 
agreed  with  the  thermocouple  temperature  to  within  10  K. 

Once  temperature  is  determined  using  either  method.  Equation  (22)  can  be  solved 
for  the  total  concentration,  Aj-,  where  the  value  of  the  partition  function  at  Tref  and  T  are 
given  by  a  piecewise,  third  degree  polynomial  fit  to  HITRAN  plots. [9,  13]  Nt  can  also  be 
determined  from  the  y-intercept  of  the  Boltzmann  plot  using  the  same  method  to  find  the 
value  of  the  partition  function  at  T. 

The  pressure  can  be  determined  using  the  Lorentzian  contribution  to  the  width  of 
the  Voigt  profile  of  the  absorption  or  by  substituting  N/V=  Nt=  P/kT  from  the  ideal  gas 
law  into  Equation  (22),  where  N  is  the  total  number  of  molecules,  V  is  the  volume  in 
cm'  ,  P  is  the  pressure,  k  is  Boltzmann’s  constant,  T  is  temperature,  and  Nj  is  the  total 
species  concentration  in  molecules-cm'  .  The  pressure  broadening  term  can  be 
determined  from  the  Igor  fitting  routine  used  to  fit  the  TDLAS  measurements  to 
determine  integrated  areas.  The  pressure  is  then  given  by: 

Avi/2=6p,  (27) 
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where  Avi/2  is  the  full  width  at  half  maximum  of  the  Lorentzian  contribution  to  the  Voigt 
line  shape  and  b  is  the  pressure  broadening  coefficient  obtained  from  HITRAN.[8] 
Pressure  determination  was  not  completed  for  the  flat  flame  burner  because  pressure 
within  the  flame  was  very  close  to  the  pressure  in  the  room.  Pressure  was  calculated  for 
the  isolator.  Preliminary  results  are  shown  in  Chapter  V. 

Tomography 

Background 

While  knowing  the  frequency  integrated  absorption  coefficient,  species  density, 
temperature,  or  pressure  along  a  given  path  is  valuable  information,  the  line  integrals 
provided  by  TDLAS  measurements  do  not  provide  any  information  about  how  these 
values  are  changing  along  the  path.  Spatial  variations  in  temperature  and  pressure  along 
the  path  can  give  important  information  about  the  flow.  Using  tomography,  interpolation 
between  the  measured  beam  paths  can  provide  the  information  necessary  to  obtain  a 
spatial  image  of  the  flow  using  a  limited  number  of  beam  paths.  [14] 

Basis  Set  Formation 

Tomographic  reconstruction  of  the  flow  characteristics  gives  much  greater  insight 
into  the  problem.  Traditional  tomography  requires  too  many  lines  of  sight  for  situations 
with  limited  optical  access.  Therefore,  two  models  were  used  to  take  advantage  of  a 
priori  information  in  the  flat  flame  burner  experiment  to  reduce  the  number  of  beam 
paths  necessary. 

The  first  model  used  was  a  two-zone  model.  This  model  assumed  radial 
symmetry  of  the  flat  flame  burner  resulting  in  a  hot  disk  surrounded  by  room-temperature 
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gas.  CFD  calculations  provided  the  second  model  used.  [15]  The  CFD  used  many  of  the 
parameters  of  the  flat  flame  experiment  to  determine  a  radial  contour  of  the  gas 
temperature  and  density.  The  CFD  models  were  provided  for  a  range  of  heights,  from  3 
mm  to  9  mm  above  the  surface  of  the  burner,  and  scaled  with  respect  to  temperature  and 
water  concentration  to  represent  various  potential  operating  conditions  of  the  system. 
Every  permutation  of  height  above  the  burner,  maximum  temperature,  and  species 
concentration  was  used  to  generate  frequency  integrated  absorption  coefficient  maps. 
The  frequency  integrated  absorption  coefficient,  a ,  is: 

S  =  (28) 

The  collection  of  these  images  covering  the  range  of  expected  operating  conditions  for 
the  flat  flame  burner  forms  what  is  called  a  training  set.  A  more  thorough  development 
of  the  process  to  create  a  frequency  integrated  absorption  map  from  CFD  calculations  is 
presented  in  Appendix  B.  Example  maps  are  presented  in  Chapter  IV.  Many  images  in 
the  training  set  provide  similar  contour  information  resulting  in  a  great  deal  of 
redundancy  in  the  training  set.  Thus,  data  reduction  techniques  can  be  applied  to 
consolidate  the  information  into  a  basis  set  which  can  then  be  used  to  represent  the 
system  under  investigation.  As  mentioned  in  the  introduction,  two  algorithms  were 
applied  to  the  CFD  data  to  obtain  the  optimal  basis  set:  vector  quantization  (VQ)  and 
principal  components  analysis  (PC  A).  The  members  of  the  basis  set  resulting  from  PC  A 
are  orthogonal  and  can,  therefore,  be  linearly  combined  to  represent  any  vector  in  the 
original  training  set.  Using  the  optimal  basis  set  and  measured  absorption  data,  a  two- 
dimensional  temperature,  density,  or  pressure  image  of  the  flow  can  be  constructed.  [3] 
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The  VQ  algorithm  was  not  as  useful  for  basis  generation  as  originally  expected. 
Instead,  VQ  was  used  primarily  as  a  clustering  and  optimization  algorithm.  It  was  used 
to  optimize  the  diameter  and  center  location  for  the  two-zone  model  of  the  flat  flame 
burner  and  the  diameter  for  the  CFD  model  of  the  flat  flame  burner.  It  was  also  used  to 
find  similar  shock  structure  in  the  time  resolved  isolator  data,  allowing  averaging  of 
similar  shock  structures  to  reduce  noise  and  allow  comparison  with  the  isolator  CFD. [16] 
Although  VQ  was  not  used  for  post-processed  tomography,  it  was  used  to  create 
the  basis  functions  for  the  near  real-time  tomography  performed  during  the  isolator 
experiment.  Due  to  the  high  degree  of  redundancy  found  in  the  isolator  CFD  data, 
images  were  presorted  based  on  the  features  contained  within  them.  The  isolator  CFD 
data  consists  of  images  with  very  distinct  regions.  Conditions  in  the  isolator  vary  from  a 
fairly  uniform  flow  to  a  shock  train  when  backpressure  is  added.  The  flow  generally 
consists  of  a  uniform  low  pressure  area  before  the  shock,  the  shock  front,  the  shock  train, 
and  a  trailing  area  of  uniform  higher  pressure  after  the  shock.  The  following  figure 
shows  an  example  of  the  CFD  results  for  the  isolator.  [3] 
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Figure  5.  CFD  model  of  frequency  integrated  absorption  coefficient  versus  position  in  the 
isolator  as  backpressure  is  increased,  blue  represents  low  and  red  represents  high 
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Note  the  position  of  the  shock  front  as  backpressure  is  increased.  While  the  image  shows 
the  entire  isolator,  in  reality  only  a  portion  of  the  isolator  is  probed  using  TDLAS.  As  the 
backpressure  was  varied,  the  shock  moved  through  the  window  where  it  was  measured  by 
TDLAS.  A  reasonable  presort,  therefore,  of  the  possible  images  consisted  of  a  uniform 
flow  to  represent  the  pre  and  post  shock  regions  and  several  shapes  to  represent  possible 
contours  of  the  shock  as  it  passed  through  the  window. 

As  an  initial  guess,  the  first  basis  set  for  the  isolator  was  obtained  by  inspection  of 
the  major  features  in  the  training  set.  VQ  then  took  this  initial  guess  as  a  basis  set, 
classified  all  of  the  training  set  images  according  to  the  initial  basis  set,  then  averaged 
over  all  of  the  images  in  each  set  to  create  a  new  basis  vector  for  each  classification.  The 
algorithm  then  proceeded  to  classify  all  of  the  images  in  the  training  set  according  to  this 
new  averaged  basis  set  and  averaged  again.  This  process  continued  until  the  relative 
change  of  the  new  basis  set  from  the  previous  basis  set  was  within  a  given  threshold. 
Once  the  training  set  was  sufficiently  classified  by  VQ,  PC  A  was  used  on  the  classified 
data  to  produce  a  basis  set  for  the  projection  of  the  TDLAS  array. [6]  The  VQ  algorithm 
can  be  found  in  Appendix  C  in  the  program  Create  VQbasis. 

The  PCA  used  for  both  the  near  real-time  tomography  as  well  as  the  post¬ 
processing  tomography  was  accomplished  by  first  finding  the  mean  of  the  training  set. 
Once  the  mean  was  determined,  the  projection  matrix  for  the  lines  of  sight  from  the  mean 
image  were  calculated  and  subtracted  from  the  data.  The  principal  components  of  the 
training  set  were  then  fit  to  the  remaining  data.  The  principal  components  of  the  training 
set  were  determined  by  finding  the  eigenvalues  and  eigenvectors  of  the  covariance  matrix 
of  the  training  set.  The  result  was  a  rank  ordered  set  of  images  where  the  highest 


25 


eigenvalue  corresponded  to  the  image  which  captured  the  highest  amount  of  variance. [5] 
The  basis  vectors  in  each  new  PCA  basis  set  are  orthogonal  to  each  other  and  can 
therefore  be  added  or  subtracted  in  linear  combinations  to  reconstruct  any  image  in  the 
training  set.  The  mean  image  and  the  first  few  principal  components  of  the  training  set 
became  the  basis  which  was  used  to  fit  the  tunable  diode  laser  data.  The  PCA  algorithm 
can  also  be  found  in  Appendix  C  in  the  Create_basisv2  program. 

Fitting  Basis  Sets  to  the  Data 

Using  PCA  to  perform  tomography  with  the  flat  flame  burner  data  allowed  a 
comparison  with  the  values  obtained  from  traditional  thermocouple  measurements.  Once 
a  basis  set  was  acquired,  the  integrated  absorptions  were  fit  to  a  linear  combination  of  the 
basis  set  using  a  least  squares  fitting  algorithm.  The  projection  matrix  for  each  basis 
image  consists  of  the  path  integral  of  each  line  of  sight  in  the  basis  image  as  shown  in 
Figure  6. [14] 


Figure  6.  Illustration  of  path  integrals  of  lines  of  sight  in  a  basis  image 


The  relation  below  transforms  the  path  to  Cartesian  coordinates  as  a  function  of  0  and 
perpendicular  distance  from  the  origin,  s: 
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p(s,0)j  =-  [ln[  dv  =  \  \a{x,y)5{s-xcos,9-ys,m.0)dxdy  =  y{ia  (29) 

L  yh{s,o,v))  II 

where  S(x)  is  the  Dirac  delta  function,  91,.  is  the  Radon  transform  of  the  /**'  projection,  and 

V  is  the  frequency  of  the  TDLAS  measurement.  [17]  Each  beam  path  can  then  be  made 
discrete  by  the  following  equation: 

Y.Vj=Pi 

i=i 

where  i  =  and  M  is  the  number  of  pixels  in  the  reconstructed  image.  In  many 

applications,  such  as  MRI  or  CAT  scans,  a  large  number  of  beam  paths  would  be 
required  to  obtain,  using  Equation  (30),  any  sort  of  spatial  resolution.  However,  using  a 
priori  information,  like  the  CFD  models,  assumptions  can  be  made  about  how  contours  in 
the  flow  will  typically  look.  CFD  models  produce  very  large  training  sets  representing 
all  of  the  possible  conditions  of  the  flow.  This  presents  a  problem  because  the  number  of 
beam  paths  required  to  make  the  fit  well  defined  increases  with  the  number  of  images 
used  to  perform  the  fit.  [17]  The  size  of  the  training  sets  can  be  reduced,  however,  using 
VQ  or  PC  A.  The  key  to  reducing  the  number  of  beam  paths  is  to  expand  the  absorption 
function  in  Equation  (30)  in  terms  of  a  smaller  set  of  basis  functions  that  span  the 
image. [3]  By  assuming  that  the  absorption  at  any  point  in  the  flow  can  be  approximated 
by  a  linear  combination  of  the  basis  set  from  the  CFD,  Equation  (29)  can  be  simplified  to: 

n  n 

a{x,y)Ka{x,y)  =  Y,qjbj(x,y)^  p^  =  iR.a  « iR.a  =  ,  (31) 

7=1  7=1 

where  qj  are  the  expansion  coefficients  for  a,  bj(x,y)  are  the  basis  functions,  and  ry  is  the 
projection  of  bj,  9? .  bj.  Equation  (29)  can  then  be  written  in  the  simplified  form  shown 
below: 
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p  =  '^q. 


(32) 


where  9^  is  the  projection  matrix  with  M  columns  and  N  rows  created  by  completing  M 
path  integrals  for  basis  functions  and  p  is  the  matrix  of  line  of  sight  measurements.  A 
least  squares  approach  is  then  used  to  solve  the  matrix  equation,  determining  the 
expansion  coefficient  vector,  q.  This  results  in  a  frequency  integrated  absorption 
coefficient  map  for  a  single  transition. [3] 

Finally,  once  the  optimal  method  for  obtaining  a  basis  set  from  the  CFD  training 
set  was  achieved  for  the  flat  flame  burner  data,  the  method  was  applied  to  the  isolator 
data.  Values  obtained  in  the  temperature  and  density  maps  were  compared  to 
thermocouple  measurements  for  the  flat  flame  burner  in  Chapter  IV.  Preliminary 
tomography  results  using  a  CFD  model  for  the  isolator  are  shown  in  Chapter  V. 
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Ill  Thermocouple  Flat  Flame  Burner  Results 


Introduction 

This  chapter  presents  the  results  of  the  thermocouple  experiment  used  to 
determine  the  temperature  profile  of  a  flat  flame  McKenna  type  burner.  The  purpose  of 
this  experiment  was  to  provide  temperature  measurements  of  the  flat  flame,  obtained 
using  traditional  thermocouple  measurements,  to  compare  with  the  temperature  results  of 
the  TDLAS  and  tomography  detailed  in  Chapter  IV.  The  first  section  gives  a  brief 
historical  background  and  explanation  of  the  theory  of  thermocouples.  Next,  the 
experimental  setup  and  procedure  are  discussed  in  detail.  Then,  the  radiative  correction, 
performed  on  the  thermocouple  data  due  to  the  high  temperatures  found  in  the  flame,  is 
discussed.  When  inside  the  flame,  the  thermocouples  became  hot  enough  to  reradiate  a 
significant  amount  of  energy,  causing  errors  in  the  thermocouple  readings  of  up  to  80  K. 
A  correction  was  performed  to  obtain  the  actual  gas  temperature  in  the  flame.  Finally, 
the  raw  and  corrected  data  are  presented  and  discussed. 

Background 

Thermocouples  were  originally  discovered  by  Thomas  Seebeck  in  1821.  Every 
conductive  metal  generates  a  potential  difference  when  it  is  subjected  to  a  thermal 
gradient.  If  one  end  is  subjected  to  a  higher  temperature  than  the  other,  the  electrons  in 
the  hot  end  have  greater  energy  and  therefore  greater  velocities.  This  means  they  are 
more  likely  to  be  found  further  from  the  hot  end.  Thus,  population  is  increased  at  the 
cold  end  resulting  in  an  electric  field  and  therefore  a  voltage  difference  between  the  two 
ends.  Thermocouples  take  advantage  of  this  fact  by  joining  two  different  types  of 
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conductors  together  at  one  end  and  measuring  the  voltage  difference  between  the  two 
metals  at  the  other  end.  [18]  Using  this  method,  however,  can  only  provide  relative 
temperature  measurements.  Therefore,  modern  voltage  meters  contain  a  calibrated 
lookup  table  to  determine  the  actual  temperature  measured. 

Thermocouple  Burner  Setup 

The  burner  used  in  this  experiment  was  a  flat  flame  McKenna  type  burner.  The 
burner  consisted  of  a  60  mm  diameter  disk  of  porous  bronze  through  which  air  and  fuel 
flows.  The  fuel  mixture  was  controlled  by  a  TYLAN  RO-28  meter  at  1.31  standard  liters 
per  minute  (SLPM)  of  Methane  and  15.6  SLPM  of  Air.  This  gave  an  equivalence  ratio, 
or  ratio  of  fuel  to  air,  of  0.8.  The  equivalence  ratio  of  0.8  means  that  the  fuel  is  lean 
compared  to  the  air,  and  therefore  allows  a  complete  combustion  model  to  be  assumed 
when  mole  fractions  are  calculated  in  Chapter  IV.  The  burner  was  cooled  by  a  Neslab 
CFT-25  cooler  with  a  MD-30  pump  flowing  water  at  274  mL/min.  Temperature, 
pressure,  and  relative  humidity  in  the  lab  were  measured  by  a  NIST  traceable  Vaisala 
PTU  (pressure,  temperature,  and  humidity)  transmitter.  Two  types  of  thermocouples 
were  used  to  measure  the  temperature  profile  of  the  flame.  Type  B  thermocouples, 
consisting  of  a  platinum-rhodium  alloy,  have  a  temperature  range  of  approximately  600 
to  2000  K  and  were  used  to  measure  temperatures  inside  the  flame.  Type  K 
thermocouples  are  made  of  chromel  and  alumel  and  have  a  usable  temperature  range  of 
approximately  70  to  700  K.  Temperatures  outside  of  the  flame  and  up  to  the  boundary  of 
the  flame  were  measured  with  type  K  thermocouples. 

Three  type  B  thermocouples  with  bead  sizes  of  0.33,  0.74,  and  1.08  mm  were 
used  to  measure  the  temperature  of  the  flame  across  the  center  of  the  burner.  The  bead 
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size  is  the  diameter  of  the  sphere  joining  the  two  thermocouple  wires  together.  The  array 
of  type  B  thermocouples  with  various  bead  sizes  enabled  the  radiative  correction 
performed  in  the  next  section.  The  thermocouples  were  mounted  on  a  translation  stage 
which  allowed  them  to  be  placed  in  the  center  of  the  burner.  The  burner  was  also 
mounted  to  a  translation  stage  allowing  the  thermocouples  to  scan  the  entire  flame 
profile.  Thin  ceramic  cylinders  were  used  to  stabilize  the  thermocouples  during  the 
experiment.  Approximately  one  quarter  inch  of  each  thermocouple  was  left  exposed  in 
the  flame  to  keep  the  heating  of  the  ceramics  from  causing  any  dielectric  effect  on  the 
thermocouple  measurements.  A  diagram  and  picture  of  the  type  B  thermocouple  array  in 
the  flame  is  shown  below. 


Figure  7.  Three  type  B  thermocouples  of  varying  bead  diameter  placed  6  mm  above  the 

surface  of  the  flat  flame  burner 

Two  type  K  thermocouples  were  used  to  measure  the  burner  profile  from  outside 
the  flame  at  room  temperature  up  to  600  K  at  the  boundary  of  the  flame.  Two  more  type 
K  thermocouples  were  used  to  measure  the  temperature  of  the  input  and  output  water 
coolant  to  the  burner.  This  allowed  a  temperature  of  combustion  to  be  calculated  as 
shown  in  the  thermocouple  results  section.  Radiative  correction  was  not  necessary  for 
the  temperatures  measured  by  the  type  K  thermocouples. 
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Radiative  Correction 


In  the  thermocouple  flat  flame  burner  experiment,  temperatures  inside  the  flame 
were  high  enough  that  the  type  B  thermocouples  reradiated  a  significant  amount  of 
energy  back  into  the  surrounding  environment  as  shown  by  the  red  glow  in  Figure  7. 
This  lost  energy  was,  therefore,  not  measured  by  the  thermocouples  resulting  in 
erroneously  low  temperature  readings.  The  amount  of  energy  lost  is  dependant  on  the 
surface  area  of  the  bead  joining  the  two  wires  of  the  thermocouple.  The  temperature 
measurement  of  the  flat  flame  burner  was  therefore  performed  using  an  array  of  three 
type  B  thermocouples  of  varying  bead  diameter.  For  an  object  in  thermal  equilibrium, 
the  heat  lost  from  radiation  is  equal  to  the  heat  gained.  In  the  case  of  the  thermocouples, 
heat  was  gained  almost  completely  through  convection  and  lost  almost  completely 
through  radiation.  Using  the  equation  for  radiation  and  Newton’s  law  of  cooling,  the 
equation  for  energy  balance  on  the  bead  is  shown  below:  [19] 

qconv=Kon.ATg-T,)  .  (33) 

KonAT^-T,)  =  6ic7,j: 

In  this  equation,  qrad  and  qconv  are  the  amounts  of  heat  lost  to  radiation  and  convection 
respectively,  e  is  the  emissivity  of  the  bead,  A  is  the  surface  area  of  the  bead,  osb  is  the 
Stefan-Boltzmann  constant,  hconv  is  the  convective  heat  transfer  coefficient  between  the 
gas  flow  and  the  bead,  and  Tb,  Tg,  and  Too  are  the  temperatures  of  the  bead,  gas  and 
surrounding  room  respectively.  The  true  temperature  of  the  gas  was  calculated  by 
rearranging  the  previous  equation  into: [19] 
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r,  =  cj(m)(r/-rj)  +  r^ 

where  C  = 

^Kc  Konv 


(34) 


then  performing  a  linear  fit  to  the  three  data  points  provided  by  the  small,  medium,  and 
large  beaded  thermocouples.  Plotting  7),  versus  Cd{m){T^  -T^^)  gives  Tg  as  the 


extrapolation  to  the  y-intercept  when  d  goes  to  zero.  In  this  equation,  N  is  the  Nusselt 
number  for  the  gas  and  bead,  ktc  is  the  gas  thermal  conductivity,  and  d(m)  is  the  diameter 
of  the  bead.  This  represents  the  perfect  thermocouple  with  a  bead  diameter  of  zero.  A 
bead  diameter  of  zero  has  no  surface  area  and  therefore  cannot  radiate. 

However,  the  situation  is  made  more  complicated  by  the  Nusselt  number’s 
dependence  on  bead  diameter,  as  shown  here:  [20] 

A  =  2  +  (0.4P^  +  0.06P^ ) 


R  = 


Av 


(35) 


where  P  is  the  Prandtl  number,  R  is  the  Reynolds  number,  pg  is  the  density  of  the  gas,  U 
is  the  velocity  of  the  gas,  and  is  the  viscosity  of  the  gas.  This  additional  dependence 
was  included  in  the  emissivity  calculation. 

At  each  location  where  all  three  thermocouples  measured  the  temperature,  data 
was  plotted  according  to  Equation  (34);  and  a  linear  fit  was  performed  to  determine  the 
actual  gas  temperature.  Figure  8  shows  an  example  of  this  process  for  three  points  near 
the  center  of  the  flame;  ignoring  the  dependence  of  the  Nusselt  number  on  bead  diameter. 
For  the  plot  shown  here,  T^a  was  calculated  by  averaging  the  temperatures  from  six  sides 
of  a  cube  where  the  bottom  of  the  cube  was  the  surface  of  the  burner  and  the  remaining 
five  sides  were  the  walls  and  ceiling  of  the  room.  The  ceiling  and  surrounding  walls 
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were  assumed  to  be  300  K  and  the  surface  of  the  burner  was  assumed  to  be  800  K,  giving 
an  average  of  383  K.  The  y-intercept  gives  Tg  as  1630  K.  As  the  emissivity  and  thermal 
conductivity  are  constant  for  the  different  bead  diameters,  their  values  had  no  effect  on 
the  gas  temperature. 


Figure  8.  Example  of  thermocouple  correction  using  Equation  (34)  with  temperatures 
from  the  three  bead  diameters  near  the  center  of  the  flat  flame  burner 


Results 

Figure  9  shows  the  temperature  profde  for  the  flat  flame  burner  before  the 
correction  was  applied.  Here  the  temperature  lost  due  to  radiation  is  readily  apparent. 
The  largest  diameter  bead  recorded  the  lowest  temperature  due  to  its  larger  surface  area 
allowing  greater  radiation  back  into  the  room.  The  smallest  diameter  bead  recorded  the 
closest  to  the  actual  temperature  because  it  is  closer  to  a  perfect  thermocouple  with  no 
surface  area  and  therefore  no  radiation  of  heat  back  into  the  room. 
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Figure  9.  Burner  Profile  before  eorreetion.  Note  as  bead  diameter  increases,  the 
heat  lost  through  radiation  increases  and  the  indicated  temperature  is  decreased 

The  three  thermocouple  diameters  were  measured  to  be  0.33  ±  0.01  mm,  0.74  ±  0.01  mm, 

and  1.08  ±  0.01  mm.  Using  the  correction  from  Equation  (34)  for  a  sample  three  points, 

the  y-intercept  in  Figure  8  gives  the  temperature  that  an  ideal  point  source  thermocouple 

would  read.  Performing  this  process  across  the  entire  profile  and  taking  the  bead 

diameter  dependence  of  the  Nusselt  number  into  account  results  in  the  following 

corrected  flat  flame  burner  profile  in  Figure  10  where  the  black  linear  fit  comes  from  the 

extrapolation  of  the  plot  of  Equation  (34)  and  the  colored  lines  are  the  thermocouples 

corrected  independently.  The  emissivity  of  the  beads  was  found  to  be  0.22  which  agrees 

well  with  previously  published  results. [21]  Note  that  the  corrected  profile  gives  a 

temperature  of  1638  ±  10  K  across  the  center  40  mm  of  the  burner. 
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Figure  10.  Corrected  burner  profile  showing  the  corrected  profile  for  each  of  the  three 
bead  sizes  and  the  fit  from  the  y-intercept  of  Equation  (34) 

Temperature  at  the  center  of  the  flame  can  also  be  obtained  using  the  adiabatic 
flame  temperature  of  the  methane  and  air  mixture.  If  the  adiabatic  temperature  of  the 
combusting  fuel  and  air  mixture  is  known,  then  the  true  flame  temperature  can  be 
calculated  by  subtracting  the  heat  lost  to  the  burner.  The  primary  source  of  heat  loss  in 
this  experiment  was  due  to  the  cooling  of  the  flat  flame  burner  with  refrigerated  water. 
The  input  and  output  temperatures  of  the  water  coolant  were  measured  with  type  K 
thermocouples.  Using  the  temperature  difference  between  the  input  and  output 
temperatures  of  the  water  coolant  along  with  the  heat  capacity  of  water  allowed  the 
calculation  of  the  amount  of  heat  lost  to  the  water.  For  the  flow  rate  of  1 .3 1  SLPM  of 
methane  at  the  equivalence  ratio  of  0.8  used  in  this  experiment,  the  adiabatic  flame 
temperature  is  1997  ±  50  K.[22]  The  measured  temperature  difference  of  water  input  to 
output  of  8.16  K,  the  measured  flow  of  the  water  of  0.27  L/min,  and  the  specific  heat  of 
water  of  75.38  J/mol-K  give  the  heat  transferred  to  the  burner  as  0.15  KJ/sec.  The 
measured  flow  of  methane  of  1.31  SLPM  and  the  heat  of  combustion  of  methane  of  802 
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KJ/mol  give  the  energy  released  by  the  methane  as  0.782  KJ/sec.  The  fraction  of  energy 
removed  by  the  water  then  is  given  by  the  ratio  of  .015  KJ/sec  to  0.782  KJ/sec;  or  0.19. 
The  adjusted  flame  temperature  is  given  as  (1-0.19)  x  1997  K  =  1618  ±  50  K  for  the 
cooled  flat  flame  burner  used  in  this  experiment.  [22] 

The  temperature  profile  of  the  flat  flame  burner  obtained  through  this  traditional 
thermocouple  measurement  experiment,  then  further  confirmed  with  a  heat  of 
combustion  calculation,  provides  a  comparison  for  the  temperatures  obtained  using 
TDLAS  over  a  single  path  and  TDLAT  with  multiple  paths  in  Chapter  IV.  Also,  the  flat 
flame,  McKenna  type  burner  was  shown  to,  in  fact,  be  flat  in  terms  of  its  temperature 
profile  inside  of  the  boundary  region.  This  lends  confidence  to  the  symmetry  arguments 
used  to  create  basis  functions  for  tomography  in  the  next  chapter.  Additionally,  the 
boundary  region  temperature  profile  was  measured  at  1  mm  resolution  using  the 
thermocouple  arrays,  which  allowed  a  good  evaluation  of  CFD  modeling  over  this  region. 
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IV  Burner  TDLAS  and  Tomography  Results 

Introduction 

The  flat  flame  burner  experiment  was  used  as  a  proof  of  concept  for  the  isolator 
experiment.  The  flow  from  the  flat  flame  burner  is  similar  to  the  flow  in  the  isolator  in 
that  neither  flow  is  uniform  along  the  line  of  sight  of  the  laser.  Also,  methane  was  used 
as  the  fuel  for  the  flat  flame  burner  and  as  the  fuel  to  vitiate  the  flow  in  the  isolator.  In 
both  cases  it  provides  a  controlled  amount  of  water  in  the  flow  to  be  observed  using 
TDLAS.  The  two  experiments  are  also  similar  in  that,  although  the  flows  are  not  uniform 
along  the  lines  of  sight,  they  provide  planar  symmetry.  This  allows  a  small  array  of 
lasers  to  provide  enough  information  to  achieve  two-dimensional  tomographic  maps  of 
the  frequency  integrated  absorption  coefficient,  temperature,  and  water  concentration. 

This  chapter  walks  through  the  entire  experiment  using  tunable  diode  lasers  to 
measure  flow  characteristics  across  the  flat  flame  burner.  The  setup  is  discussed  first, 
including  an  in-depth  look  at  the  data  acquisition  system  from  the  function  generators 
controlling  the  lasers  to  the  detection  and  conversion  of  the  analog  laser  intensity  to 
digital  measurements.  Next,  the  logic  behind  the  wavelength  range  chosen  for  each  laser 
is  presented.  The  third  section  discusses  the  process  of  converting  raw  intensity  data  to 
absorption  spectra  and  integrated  absorption.  The  theory  behind  the  two-zone  analytic 
model  as  well  as  the  resulting  temperature  and  density  determination  are  presented  next. 
Finally,  the  tomography  section  includes  the  application  of  the  same  two-zone  model  as 
well  as  the  CFD  model  of  the  flat  flame  burner  with  a  final  comparison  of  all  of  the 
methods  to  the  thermocouple  experiment  results  from  Chapter  III. 
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TDLAS  Burner  Setup 


TDLAS  was  performed  on  the  same  flat  flame  McKenna  type  burner  running  at 
the  same  parameters  used  in  the  thermocouple  experiment  presented  in  Chapter  III.  The 
cylindrical  flat  flame  burner  provided  a  flow  with  elevated  water  concentration  and 
temperature  in  a  simplified  geometry  to  be  observed  by  a  laser  array.  Due  to  the  fact  that 
the  lines  of  sight  pass  through  room  temperature  air  on  either  side  of  the  flame,  the  flow 
is  not  uniform  along  the  path  of  each  laser.  Figure  1 1  shows  the  schematic  of  the  entire 
setup. 


Figure  11.  Schematic  of  optical  arrangement  for  flat-flame  burner  and  isolator 

experiments 


Note  that  the  burner  experiment  used  only  one  of  the  tunable  diode  laser  arrays  and  only 
seven  beam  paths.  A  National  Instruments  (NI)  PXI-8105  computer  ran  the  experimental 
setup  for  the  flat  flame  burner  and  isolator.  It  sent  a  control  signal  to  three  NI  PXI-5402 
14-bit,  100  MS/s  function  generators  to  provide  the  current  modulation  signal  to  an  ILX 
driver.  The  ILX  driver  provided  thermo-electric  cooling  and  current  to  the  three  lasers. 
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The  three  lasers,  operating  at  nominal  wavelengths  of  1391  nm,  1392  nm,  and  1397  nm, 
were  combined  into  a  single  fiber  by  a  4x1  Newport  combiner.  The  combined  signal 
then  continued  to  a  1x2  Newport  splitter  which  sent  10%  of  the  signal  to  a  50/50  1x2 
Newport  splitter  and  the  remaining  90%  of  the  signal  to  another  50/50  splitter.  The  50/50 
splitter  receiving  10%  of  the  signal  directed  one  signal  to  a  SiOi  etalon  for  a  wavelength 
measurement  while  the  other  was  used  as  a  reference  beam.  The  50/50  splitter  receiving 
90%  of  the  signal  sent  each  signal  to  a  separate  1x8  splitter.  The  two  1x8  splitters 
provided  up  to  16  beam  paths,  each  containing  three  wavelength  ranges  for  a  total  of  48 
water  absorption  spectra.  The  fibers  were  coupled  to  collimators  which  direct  collimated, 
1  mm  diameter  beams  through  the  flame  to  be  recorded  on  the  other  side  by  an  in-house 
designed  InGaAs  detector  array  complete  with  individual  preamplifiers.  An  overhead 
image  of  the  laser  array  and  detector  setup  for  the  flat  flame  burner  experiment  is  shown 
in  Figure  12. [23]  The  recorded  analog  spectra  are  sent  from  the  board  to  a  NI 
analog/digital  converter  before  they  are  sent  back  to  the  NI  PXI-8105  computer  to  be 
recorded  into  memory.  [3]  The  entire  setup  enables  near  real-time  data  capture  and 
tomography. 
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Figure  12.  Overhead  view  of  flat  flame  burner  experiment  with  lines  of  sight  overlaid 

Selection  of  Observed  Transitions 

Three  lasers  were  used  for  the  flat  flame  burner  experiment.  Rough  tuning  for 
each  of  the  three  lasers  was  performed  using  a  thermal  electric  cooler  to  obtain  nominal 
wavelengths  of  1391,  1392,  and  1397  nm.  This  allowed  different  wavelength  ranges,  and 
therefore  different  transitions  of  the  water  molecule,  to  be  observed  by  each  laser.  Table 
1  and  the  corresponding  graphs  in  Figure  13  show  the  identified  water  transitions  in  the 
wavelength  regions  used. 
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Table  1.  The  12  transitions  used  for  experiment  are  shown  in  order  of  increasing 


wavelength 


X  (nm) 

Vo  (cm'‘) 

S  at  300  K  (cm'- 
molecule'-cm'^)x  10'^^ 

Eo  (cm"') 

v' 

fKa’Kc’ 

Label 

1391.67 

7185.60 

83.35 

1045.06 

1  0  1 

6  6  0 

6 

6  1 

a 

1391.71 

7185.40 

1.131 

1474.98 

1  0  1 

9  5  4 

9 

5  5 

b 

1392.53 

7181.16 

1488 

136.76 

1  0  1 

2  0  2 

3 

0  3 

c 

1392.55 

7181.10 

1211> 

1216.20 

120 

7  7  0 

7 

6  1 

d 

1392.64 

7180.61 

3.277 

1477.30 

1  0  1 

9  5  5 

9 

5  4 

e 

1392.68 

7180.40 

55.77 

224.84 

200 

3  2  1 

4 

1  4 

f 

1392.81 

7179.75 

24.38 

1216.20 

1  0  1 

7  6  2 

7 

6  1 

g 

1396.37 

7161.41 

1168 

224.84 

1  0  1 

3  1  3 

4 

1  4 

h 

1396.42 

7161.19 

.7382 

1446.13 

02  1 

11  3  9 

10 

3  8 

i 

1396.49 

7160.81 

.5983 

1874.97 

1  0  1 

10  6  4 

10 

6  5 

j 

1396.60 

7160.24 

421.2 

222.05 

1  0  1 

3  0  3 

4 

0  4 

k 

1396.61 

7160.19 

121.9 

173.37 

200 

3  0  3 

3 

1  2 

1 

corresponding  to  Table  1 


The  selection  of  these  wavelength  ranges  was  based  on  the  specific  peaks  captured  by 
each  laser.  Peaks  were  chosen  first  based  on  their  intensity,  where  higher  intensity  allows 
them  to  be  integrated  more  accurately.  Second,  peaks  were  chosen  to  have  a  range  of 
term  energies.  This  was  important  to  temperature  determination  because  the  x-axis  of  the 
Boltzmann  plots  is  term  energy.  Thus,  to  provide  confidence  in  the  linear  fit  to  the 
Boltzmann  plots,  a  range  of  initial  state  energies  was  required. 

Igor  Pro  was  used  to  perform  integrations  of  all  12  transitions.  [7]  The  fitting 
routine  implements  the  Humlicek  algorithm  to  obtain  a  constrained,  multi-peak  fit.  [24] 
This  method  provided  six  data  points  for  the  Boltzmann  plots  resulting  in  high  accuracy 
temperature  calculations.  These  results  are  shown  in  the  section  covering  the  two-zone 
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analytic  model.  Integrating  eaeh  spectrum  individually  becomes  time  eonsuming  when 
many  beam  paths  are  considered.  For  the  number  of  beam  paths  necessary  to  allow 
accurate  tomography,  an  automated  integration  routine  was  developed  in  Igor.  The  three 
peaks  with  the  strongest  transitions,  and  therefore  the  lowest  error  for  the  automated 
integration,  were  the  225  cm'',  the  1045  cm'',  and  the  1216  cm  '  transition.  These  peaks 
were  used  for  the  tomography  results  shown  later  in  this  chapter.  Using  fewer  peaks 
results  in  a  less  aceurate  temperature  profile,  however,  so  the  development  of  a  more 
rigorous  data  analysis  routine  that  could  fit  more  peaks  automatieally  would  increase  the 
accuraey  of  the  temperature  and  density  results. 


From  Raw  Data  to  Integrated  Absorption 

The  light  intensity  observed  at  the  deteetor  is  eonverted  into  an  absorption 
speetrum,  using  Equation  (22),  by  taking  the  natural  logarithm  of  the  deteeted  beam 
divided  by  the  reference  beam.  Figure  14  shows  one  millisecond  of  raw  data  and  the 
eorresponding  absorption  speetra  for  each  laser. 


Figure  14.  A  single  period  of  the  raw  data  from  the  three  multiplexed  lasers  with  the 

corresponding  absorption  speetra 
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A  multi-peak,  baseline  Voigt  fit  was  performed  using  Igor  to  fit  the  speetra.  The 
coeffieients  of  the  fit  were  then  used  to  determine  the  areas  of  eaeh  peak  in  the  spectrum. 
The  results  are  the  frequency  integrated  peaks,  in  units  of  wavenumbers,  which  are  equal 
to  the  right  hand  side  of  Equation  (22).  The  DAQ  software  developed  to  run  the  flat 
flame  experiments  can  perform  integrations  over  all  of  the  beam  paths,  seven  for  the  flat 
flame  burner,  of  a  single  selected  peak  in  near  real-time.  The  seven  lines  of  sight  are 
shown  below  with  the  burner  in  the  center  position. 


Figure  15.  Schematic  showing  the  seven  lines  of  sight  used  for  the  flat  flame  burner 

experiment 

The  resulting  plots  in  Figure  16  show  the  frequency  integrated  absorption  for  the  1216 
cm''  transition  for  the  seven  paths  of  the  flat  flame  burner  experiment  with  the  burner  in 
three  positions.  As  the  burner  is  moved  left  and  right,  the  profile  changes,  allowing  flame 
position  tracking  to  within  millimeter  precision.  With  only  seven  beam  paths,  the 
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information  necessary  to  determine  positions  and  diameter  to  within  millimeters  is 
available.  Similar  plots  will  be  shown  later  for  the  isolator  results  in  Chapter  VI. 


Figure  16.  Integrated  absorption  of  the  1216  cm  '  transition  for  each  of  the  seven  lines  of 
sight  with  the  flame  9  mm  left  of  center,  in  the  center  position,  and  9  mm  right  of  center 

Two  major  point  of  interest  are  apparent  in  the  graphs  in  Figure  16:  flame  position 
and  line  of  sight  path  length.  First,  the  1216  cm  '  transition  has  a  smaller  initial  state 
population  outside  of  the  flames.  This  is  shown  in  the  center  graph  where  both  lines  of 
sight,  zero  and  six,  are  completely  outside  of  the  flame.  As  the  flame  is  shifted  left  or 
right,  the  lines  of  sight  show  increased  frequency  integrated  absorption  {al)  in  units  of 
cm  '  when  the  flame  moves  into  their  respective  paths.  Differences  in  path  length  can 
also  be  extracted  from  these  graphs.  Using  the  center  position,  in  the  middle  graph  in 
Figure  16  and  in  the  schematic  of  Figure  15,  paths  one  through  five  are  shown  to  pass 
through  approximately  the  same  length  of  hot  gas,  so  the  contribution  to  the  intensity 
from  the  flame  is  about  the  same.  Paths  one,  three,  and  five,  however,  show  higher 
intensities  than  paths  two  and  four.  This  is  due  to  the  longer  cold  path  length  of  paths 
one,  three,  and  five.  Even  though  the  contribution  to  the  intensity  for  the  1216  cm  ' 
transition  from  the  cold  region  is  not  as  great  as  the  contribution  from  the  flame  region, 
the  extended  path  lengths  of  paths  one,  three,  and  five  provide  enough  added  absorption 
that  the  intensity  is  noticeably  increased.  The  information  available  in  the  graphs  of 
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Figure  16,  for  all  three  transitions,  enables  accurate  determination  of  the  temperature  both 
inside  and  outside  of  the  flame,  as  shown  in  the  next  section. 


Analytic  Two-Zone  Model 

Before  any  peaks  are  integrated,  information  can  be  extracted  about  the  burner 
using  a  two-zone  analytic  model  by  breaking  spectra  into  their  constituent  hot  and  cold 
contributions.  Using  the  assumption  that  the  burner  experiment  consists  of  a  hot  circular 
region  surrounded  by  room  temperature  air,  a  geometric  model  describing  the  system  can 
be  applied.  First,  consider  a  beam  passing  through  a  region  of  room  temperature  gas 
followed  by  a  region  of  hot  gas  before  hitting  the  detector  as  shown  below  in  Figure  17. 
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Figure  17.  Diagram  showing  absorption  of  light  through  a  cold  and  hot  region  of  gas 
The  intensity  of  the  incident  light  after  passing  through  the  two  regions  of  gas  is 
governed  by  Beer's  law: 


j  -  T 

il  — 


-S\T^)NrJ2g2i^) 


I2  — 

Substituting  in  for  /;  for  the  second  region  gives: 


(36) 


1 2  I 


-S'(71)Vr,;igi(v)  -S'(T2)Nr^l2g2(v)  _  j  -(S'(7’i)Vj.,;,g,(v)+S’(7’2)Vr,/2g2(v)) 
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Dividing  by  Iq  and  taking  the  natural  logarithm  gives: 

=  S\T,  kg,  (v)  +  S\T2  )Nrj2g2  (v)  • 


In 


J 


(38) 
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Now,  considering  the  left  hand  side  of  Equation  (38)  as  the  total  spectrum  over  the  entire 
path,  ln(Io/l2)  can  be  labeled  as  At.  The  two  terms  on  the  right  hand  side  of  Equation  (38) 
are  the  contribution  to  the  spectrum  from  the  cold  region  and  the  hot  region  respectively, 
so  S'  (Ti)NTiligi(v)  can  be  labeled  as  Ac,  and  S'  (T2)NT2l2g2(v)  can  be  labeled  as  Ah. 
Dividing  Ac  by  the  cold  path  length,  //,  and  Ah  by  the  hot  path  length,  I2,  then  gives  the 
cold  and  hot  spectra  per  unit  length.  Using  the  beam  geometry  shown  in  Figure  18,  the 
outside  beam  paths  can  be  averaged  and  divided  by  the  path  length  to  give  Ac/l}.  To  find 
Ah/I 2,  however,  the  diameter  must  first  be  determined.  Once  these  spectra  per  unit  length 
are  determined,  they  can  then  be  used  as  reference  spectra  and  multiplied  by  the  lengths 
of  the  hot  and  cold  spectra  of  any  beam  path  to  reconstruct  the  total  spectrum  as: 

(39) 

Lc  J-J 

where  Lc  is  the  path  length  of  the  cold  reference  spectra,  previously  called  h,  D  is  the 
diameter  of  the  flame,  L  is  the  path  length  of  the  line  of  sight  under  investigation,  and  d  is 
the  path  length  of  the  line  of  sight  through  the  flame,  previously  called  h.  The  change  of 
variables  allows  the  diameter  to  be  calculated  using  a  cold  reference  path,  with  no  hot 
contribution,  and  two  paths  with  different  lengths  of  hot  contribution.  [14,  17]  The  setup 
described  here  is  shown  in  Figure  18. 
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Figure  18.  Diagram  of  burner  with  line  of  sight  overlay 
The  diameter  is  determined  by  solving  for  Ah  of  two  paths  with  different  path 
lengths  through  the  flame,  d,  and  setting  them  equal  to  each  other.  The  result  after  some 
algebra  is: 


(40) 

d  Aj—  A^L 

Rearranging  Equation  (40)  enables  the  solver  in  Microsoft  Excel  to  determine  the  ratio  of 
d'/d  which  minimizes  the  following  equation: 


{A^'-A,L)-^{A,-A,L)=Q.  (41) 

d 

The  diameter  can  then  be  determined  using  geometry,  where 
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(42) 


d'  _2yl(D/2f-x'^ 
d  2-s]iD/2y-x^  ' 

Once  D  is  known,  the  cold  contribution  to  the  center  path  length  can  be  subtracted 
off,  leaving  a  reference  hot  spectrum.  Ah.  With  the  reference  spectra  Ac  and  Ah  and  their 
corresponding  lengths  Lc  and  D  determined,  all  of  the  lines  of  sight  can  be  split  into  their 
hot  and  cold  spectrum  contributions.  This  allows  for  more  accurate  integration  of  the 
absorption  spectra  which  enables  very  accurate  path  length,  and  temperature 
determination. 

Because  the  initial  level  populations  for  each  transition  are  dependent  on 
temperature,  they  give  different  diameters  of  the  flame.  For  example,  due  to  its  lower 
term  energy,  the  initial  state  of  the  1045  cm"'  transition  (in  the  1391  nm  spectrum)  begins 
to  fill  thermally  before  the  initial  state  of  the  1216  cm"*  transition  (in  the  1392  spectrum). 
Thus,  as  the  line  of  sight  moves  closer  to  the  edge  of  the  burner  and  temperature 
increases,  the  frequency  integrated  absorption  coefficient  for  the  1045  cm  *  transition 
becomes  larger  before  the  frequency  integrated  absorption  coefficient  for  the  1216  cm"* 
transition.  The  result  is  a  larger  diameter  determination  when  using  the  1045  cm"* 
transition.  The  diameter  results  for  each  spectrum  are  shown  in  Table  2. 


Table  2.  Results  of  diameter  determination  using  two-zone  analytical  method 


TDL  Center  Wavelength  (nm): 

1391 

1392 

1397 

Corresponding  Transition  (cm"*): 

1045 

1216 

224 

Diameter  (mm): 

66.9 

62.0 

56.9 

With  separate  hot  and  cold  spectra,  the  integrations  of  the  various  lines  can  be 
performed  very  accurately.  Using  the  separated  hot  and  cold  spectra,  where  the  diameter 
of  the  flame  was  determined  with  the  1397  nm  spectrum,  Boltzmann  plots  were  created 
for  the  hot  and  cold  region,  as  shown  in  Figure  19. 
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1391  nm  1392  nm  1397  nm 


Figure  19.  Separated  spectra  and  Boltzmann  plots  resulting  from  peak  integrations 
The  resultant  temperatures  and  calculated  mole  fractions  compare  extremely  well 
with  the  thermocouple  measurements  as  shown  in  Table  3.  In  the  table,  the  40  and  58 
mm  averages  refer  to  the  length  across  the  flame  used  to  calculate  an  average 
temperature.  For  example,  the  40  mm  average  uses  the  center  40  mm  of  the  temperature 
profile  shown  in  Figure  10,  from  approximately  14  mm  to  54  mm.  Note  the  mole 
fraction  for  the  40  mm  average  of  the  hot  region,  0.1547,  comes  from  a  model  assuming 
complete  combustion.  [22]  This  number  was  not  available  for  the  58  mm  average  of  the 
thermocouple  data  which  includes  some  of  the  boundary  layer. 


Table  3.  TDLAS  Temperature  and  mole  fraction  results  and  comparison  to  thermocouple 
measurements  and  calculated  mole  fractions. 


Measurement 

Temperature  (K) 

Mole  Fraction 

TDLAS  Hot 

1628  ±  20 

0.157+0.008 

TC  Hot  (58  mm  average) 

1583  ± 124 

N/A 

TC  Hot  (40  mm  average) 

1638  ± 10 

0.1547 

TDLAS  Cold 

300  ±1 

0.0071+0.0004 

TC  Cold  (at  15  mm  from  flame) 

302  ±2 

0.0061±0.0004 
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The  results  from  the  section  are  obtained  using  a  completely  analytical  approach 
to  perform  tomography.  In  a  flow  where  conditions  are  changing,  like  that  of  a  scramjet, 
many  more  datasets  would  need  to  be  acquired,  resulting  in  more  data  and  necessitating  a 
computational  approach  to  perform  tomography.  This  same  two-zone  model  can  still  be 
used  in  a  computational  approach,  as  shown  in  the  next  section. 

Computational  Tomography 

Tomography  for  the  flat  flame  burner  was  performed  using  two  models  for  basis 
set  generation.  The  first  method  uses  the  same  basis  set  as  the  analytical  two-zone  model 
where  a  radially  symmetric  hot  center  region  is  surrounded  by  cooler,  room-temperature 
gas.  It  uses  beam  paths  outside  of  the  flame  to  define  the  cold  region,  and  then  performs 
a  fit  to  capture  the  flame  region.  The  second  model  uses  a  radial  temperature  and  density 
profile  from  CFD.  This  model  provides  contour  information  in  the  boundary  region, 
which  the  two-zone  method  cannot  capture. 

Two-zone  Model 

The  two-zone  model  is  a  very  basic  example  of  the  use  of  a  basis  created  from  a 
priori  knowledge  to  reduce  the  number  of  lines  of  sight  necessary  to  perform 
tomography.  It  uses  the  same  symmetry  arguments  as  the  two-zone  analytic  model  and 
therefore  provides  an  excellent  comparison.  The  three  transitions  with  the  smallest 
integration  error  were  used  for  the  tomography  portion  of  the  experiment:  the  225  cm'*, 
1045  cm'\  and  1216  cm'*  transition.  For  each  transition,  an  absorption  coefficient  map 
was  calculated  by  first  constraining  the  cold  region  using  only  the  beam  paths  outside  of 
the  flame,  then  fitting  the  hot  region  separately.  See  the  Tophatfixedcold  program  in 
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Appendix  C  for  the  actual  code.  To  fit  the  data  to  the  hot  region,  the  projection  matrix 
for  the  beams  in  the  cold  region  was  used  to  determine  the  amount  of  integrated 
absorption  that  contributed  to  each  line  of  sight  due  to  the  cold  region.  These  path 
integrals  were  then  subtracted  from  the  integrated  absorption  given  by  the  detectors.  The 
remaining  absorption  was  then  due  completely  to  the  hot  region  which  was  fit 
independently.  The  independent  hot  region  fit  was  done  by  making  an  absorption  map 
where  the  hot  region  was  set  to  one  and  the  cold  region  was  set  to  zero  as  shown  in 
Figure  20.  The  adjusted  data  were  fit  to  the  hot  basis  map  by  scaling  it  with  a  coefficient. 
The  fitted  integrated  absorption  map  was  the  result  of  the  fixed  cold  basis  map  added  to 
the  scaled  hot  basis  map. 


Cold  Basis  Hot  Basis  Fitted 


Figure  20.  Basis  images  and  resultant  fit  absorption  map 
The  automation  of  the  fitting  process  lends  itself  to  an  optimization  where  the 
flame  diameter  and  the  center  location  of  the  basis  set  are  varied.  The  residual  of  the  fit, 
as  these  parameters  are  varied,  shows  how  well  the  basis  model  matches  the  data 
collected  at  the  detectors  using  the  following  equation:  [3,  16] 

#LOS 

Residual  =  ^  abs\a^L-  -  C  *  /? J ,  (43) 

1=1 


where  is  the  frequency  integrated  absorption  coefficient  from  the  i*  beam  path,  Li  is 
the  length  of  the  beam  path,  C  is  the  fit  coefficient  of  the  basis  image,  and  pi  is  the  line 
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integral  along  z*  beam  path  in  the  basis  image.  Using  this  logic,  the  minimum  values  of 
the  residual  maps  were  used  to  determine  where  the  best  fit  occurred.  The  diameter 
determined  using  the  optimization  technique  matches  up  well  to  the  measured  flame 
diameter,  as  well  as  the  diameter  given  by  the  two-zone  analytic  model. 


Two-Zone  Optimization 


Diameter  (cm)  8  „  in  12 

X  Position  (cm) 

Figure  21.  Two-zone  optimization  showing  residual  as  a  function  of  x  position  of  the 
center  of  the  burner  and  the  flame  diameter  are  varied 

The  plot  in  Figure  21  used  63  lines  of  sight,  so  both  diameter  and  x  position  are 
well  defined  resulting  in  a  large  curvature  at  the  minimum.  Table  4  shows  the  measured 
diameter  of  the  porous  surface  of  the  flat  flame  burner  as  well  as  the  analytic  two-zone 
model  solution  for  the  diameter  using  each  of  the  three  transitions.  Table  5  shows  the 
results  of  all  of  the  different  combinations  of  transitions  and  number  of  beams  used.  All 
of  the  residual  based  center  x  position  results  agree  with  the  measured  center  of  10  cm  to 
within  0.3  cm.  The  diameter  results  agree  with  the  physical  diameter  of  the  burner  to 
within  0.5  cm.  Also  note  the  correlation  between  the  analytic  method  diameters  and  the 
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tomography  residual  diameters  for  each  transition.  Due  to  their  sensitivities  at  different 
temperatures,  each  of  the  transitions  begins  to  pick  up  the  boundary  layer  at  different 
points.  Both  methods  show  the  225  cm'*  transition  obtaining  the  smallest  diameter  of  5.7 
±  0.2  cm.  The  next  largest  diameter  determination  comes  from  the  1216  cm”'  transition, 
with  both  methods  giving  6.1  ±  0.3  cm.  The  1045  cm''  transition  gave  the  largest 
diameter  of  6.4  ±  0.2  cm. 


Table  4.  Measured  and  calculated  flame  diameters  from  physical  apparatus  and  two-zone 

analytic  model 


Transition  (cm  '): 

225 

1045 

1216 

Physical  Diameter  of  burner  (cm): 

_ 6,0 _ ^1 

Anal54ic  Diameter  (cm): 

5.7 

6.7 

6:2  1 

Table  5.  Results  of  two-zone  optimization  of  the  x  position  of  the  center  of  the  burner  and 

the  diameter  of  the  flame 


Transition  (cm”') 

Beams  Used 

X  Position  (cm) 

Diameter  (cm) 

7 

10.3 

5.8 

225 

63 

10.2 

5.5 

119 

10.2 

5.5 

7 

10.0 

6.5 

1045 

63 

10.0 

6.5 

119 

10.0 

6.5 

7 

10.1 

5.8 

1216 

63 

10.1 

6.3 

119 

10.1 

6.3 

Due  to  the  large  number  of  beam  paths  crossing  the  x  direction,  as  shown  in 
Figure  18,  there  is  enough  information  to  provide  a  good  estimate  of  the  center  position. 
In  the  y  direction,  however,  information  is  provided  only  by  the  diagonal  beam  paths. 
The  surface  plot  of  the  two-zone  residual  map  as  both  the  x  and  y  center  positions  are 
varied  as  shown  in  Figure  22. 
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TwOnZone  Optimization 


X  Position  (cm)  Y  Position  (cm) 

Figure  22.  Two-zone  optimization  of  the  x  and  y  position  of  the  center  of  the  burner 

using  tomography 


While  the  curvature  at  the  minimum  of  the  residual  map  in  the  y  direction  is  not 
as  large  as  the  x  direction,  a  clear  minimum  is  still  discerned.  Table  6  shows  the  results 
of  each  of  the  transitions  with  three  different  numbers  of  path  lengths.  Even  with  the 
limited  information  in  the  y  direction,  the  center  y  position  is  located  to  within  1  cm  of 
the  measured  position  of  16.7  cm. 


Table  6.  Results  of  tomography  optimization  to  find  the  center  position  of  the  flame 


Transition  (cm'*) 

Beams  Used 

X  Position  (cm) 

Y  Position  (cm) 

Measured 

- 

10.0 

16.7 

7 

10.4 

17.1 

225 

63 

10.2 

16.3 

119 

10.2 

16.4 

7 

10.0 

16.7 

1045 

63 

10.1 

17.4 

119 

10.0 

17.3 

7 

10.1 

16.4 

1216 

63 

10.1 

17.3 

119 

10.1 

17.1 
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Once  a  fitted  integrated  absorption  map  was  completed  for  each  of  the  three 
transitions,  a  Boltzmann  plot  was  then  constructed  pixel  by  pixel.  From  the  temperature 
plot,  a  water  density  was  calculated.  The  temperature  and  density  surface  plots  using  the 
optimized  diameter  from  the  1216  cm"'  line  of  6.3  cm  are  shown  here  in  Figure  23.  A 
comparison  to  the  thermocouple  results  and  CFD  tomography  results  is  shown  at  the  end 
of  the  next  section. 
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Figure  23.  Temperature  and  density  maps  using  tomography  with  the  two-zone  model 

basis  images 


CFD  Model 

The  CFD  calculations  provide  the  temperature  and  density  as  a  function  of 
distance  from  the  center  of  the  flame  at  varying  heights  above  the  burner  surface.  Radial 
symmetry  was  assumed  to  create  temperature  and  density  maps.  The  radius  was  scaled 
using  an  optimization  technique  similar  to  that  used  in  the  two-zone  model.  See  the 
program  Bestratio  in  Appendix  C  for  the  actual  algorithm.  The  radial  plots  of  the  CFD 
calculations  are  shown  in  Figure  24.  Again,  note  the  difference  in  sensitivity  of  the  three 
transitions  to  the  boundary  layer.  These  differences  result  in  the  different  diameter 
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Freq.  Int.  absorption  coef.  (cm'  )x10’  Density  (Molecules/cm  )x10 


determinations  of  both  the  two-zone  analytical  model  as  well  as  the  two-zone 
optimization. 


Radial  Distance  (cm) 


Figure  24.  Radial  temperature  and  density  from  CFD  calculations  and  the  resulting  radial 
absorptions  for  the  three  transitions  used  for  tomography 
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Temperature  (K) 


The  results  of  the  radius  optimization  are  shown  in  Figure  25.  A  fit  residual  was 
calculated  for  each  of  the  three  transitions  over  a  range  of  scaled  radii  of  the  CFD  from 
0.75  to  1.25  in  increments  of  0.01.  The  range  over  which  a  good  fit  to  the  basis  functions 
can  be  obtained  is  different  for  each  transition.  Also,  the  residuals  are  fairly  flat  for  a 
range  of  values  in  which  a  good  fit  can  be  achieved  by  scaling  the  diameter  and 
temperature  of  the  flame.  Essentially,  a  basis  image  with  a  smaller  flame  diameter  can 
still  maintain  the  same  line  integral  to  fit  the  data  as  long  as  the  temperature  of  the  flame 
is  increased  at  the  same  time.  However,  once  the  diameter  becomes  so  large  that  it 
pushes  into  the  path  of  a  line  of  sight  that  is  showing  no  flame  contribution,  the  residual 
rises  dramatically.  Conversely,  if  the  diameter  becomes  so  small  that  a  line  of  sight 
which  does  show  flame  contribution  is  left  outside  of  the  flame,  the  residual  quickly 
increases  again.  To  overcome  these  difficulties,  a  weighted  average  of  the  best  fit  radius 
was  performed  by  adding  the  inverse  of  the  three  individual  residual  plots.  This  method 
weights  the  residual  plots  with  the  lowest  residual,  and  therefore  best  fit,  the  highest.  The 
result  is  shown  by  the  black  line  trace  in  Figure  25.  The  maximum  of  the  weighted  sum 
of  the  inverse  residuals  is  achieved  at  a  value  of  0.9.  This  indicated  that  the  x-axis  of  the 
radial  temperature  and  density  profiles  given  by  the  CFD  calculations  needed  to  be  scaled 
by  a  factor  of  0.9  to  fit  the  TDLAS  data. 
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CFD  Radius  Optimization 


Weighted  Combined  1 /Residual  -  225  cm'’  Residual 

-  1045  cm"’  Residual  - 1216  cm"’  Residual 


Figure  25.  Residuals  and  combined  weighted  inverse  of  residuals  showing  the 
optimization  of  the  radius  of  the  CFD  training  set 

The  training  set  was  constructed  using  the  CFD  calculations  over  a  range  of 
heights  from  3  mm  to  9  mm  above  the  surface  of  the  burner.  Each  of  the  heights  was 
then  scaled  to  a  maximum  temperature  and  a  peak  density  to  minimum  density  ratio.  The 
resulting  training  set  consisted  of  175  integrated  absorption  maps  which  capture  a  wide 
range  of  possible  conditions  of  the  burner.  Using  all  175  snapshots  to  fit  the  data  for  each 
line,  however,  would  make  the  fit  underdetermined  unless  more  lines  of  sight  were 
added.  To  achieve  a  well-posed  problem,  principal  component  analysis  was  used  to 
reduce  the  dimensionality  of  the  training  set  and  thereby  reduce  the  number  of  lines  of 
sight  necessary. 

The  first  step  of  the  PCA  algorithm  fit  the  data  to  the  mean  of  the  training  set. 
Then,  the  projection  matrix  from  the  mean  image  was  subtracted  from  the  data  like  the 
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two-zone  model  where  the  fixed  cold  contribution  was  subtracted  from  the  data.  The 


remaining  integrated  absorption  coefficient  for  each  path  was  fit  to  a  linear  combination 
of  a  justified  number  of  the  first  principal  components. 

The  justification  for  the  number  of  principal  components  used  was  provided  by 
the  representation  error  and  the  condition  number.  The  representation  error  is  a 
percentage  error  that  describes  how  well,  on  average,  each  image  in  the  training  set  can 
be  constructed  as  a  function  of  the  number  of  principal  components  kept  to  do  the 
reconstruction.  See  Reperr2  in  Appendix  C  for  the  algorithm.  For  complex  distributions, 
keeping  more  images  in  the  basis  set  allows  better  representation  of  every  image  in  the 
training  set.  The  condition  number  can  be  thought  of  as  an  error  multiplier.  The  error  in 
the  measured  data  will  be  greater  in  the  fit  coefficients  by  a  factor  of  approximately  the 
condition  number.  Condition  numbers  around  10  or  less  were  used  to  define  a  well  posed 
mathematical  problem  for  this  experiment. 


Table  7.  Condition  number  under  varying  conditions 


Eo  (cm"') 

7  Beampaths 

63  Beampaths 

1 19  Beampaths 

3  BASIS 

225 

148 

100 

40 

VECTORS 

1045 

528 

26 

10 

1216 

324 

13 

7 

2  BASIS 

225 

31 

16 

12 

VECTORS 

1045 

9 

4 

4 

1216 

4 

4 

5 

1  BASIS 

225 

1 

1 

1 

VECTOR 

1045 

1 

1 

1 

1216 

1 

1 

1 

Table  8.  Average  representation  percent  error  for  the  three  transitions  using  varying 
_ numbers  of  principal  components _ 


225  cm‘* 

1045  cm"' 

1216  cm  * 

3  BASIS  VECTORS 

.03  % 

.13  % 

.20  % 

2  BASIS  VECTORS 

.18% 

.25  % 

.34  % 

1  BASIS  VECTORS 

.46  % 

8.34  % 

17.20% 
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As  shown  in  Table  7,  using  more  lines  of  sight  lowers  the  condition  number 
because  additional  information  is  available  to  constrain  the  fit.  Table  8  shows  that  the 
representation  error  is  under  one  percent  when  two  principal  components  are  retained. 
Using  the  information  of  Table  7  and  Table  8  together,  keeping  two  principal  components 
with  1 19  beam  paths  was  chosen  as  the  optimum  solution.  Figure  26  shows  the  mean  and 
first  two  principal  components  for  the  1216  cm  '  transition.  The  linear  combination  of 
these  three  images  results  in  the  fit  absorption  map. 


Mean  Image 

-3 


0  0 


y  position  (cm) 


X  position  (cm) 


Figure  26. a.  The  mean  image  of  the  training  set  for  the  1216  cm  '  transition 
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First  Principal  Component 
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Figure  26.b.  The  first  principal  component  for  the  1216  cm  *  transition 


Second  Principal  Component 
-4 


X  10 


y  position  (cm)  ^  (cm) 

Figure  26.c.  The  second  principal  component  for  the  1216  cm  *  transition 
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Frequency  Integrated  Absorption  Coefficient  Map 
Kin'' 


0  0 

y  position  (cm)  position  (cm) 

Figure  26.d.  Fit  absorption  map  for  the  1216  cm  '  transition 

Figure  26.  Fit  coefficients  for  the  mean  image  and  first  two  principal  components  are 
multiplied  by  their  corresponding  basis  image  then  added  in  linear  combination  to 
produce  the  fitted  frequency  integrated  absorption  coefficient  map 

Like  the  two-zone  model,  temperatures  were  calculated  pixel  by  pixel  using  the 
three  fitted  integrated  absorption  maps  while  density  was  calculated  from  the  temperature 
map  using  Equation  (18).  Figure  27  shows  the  absorption  maps  using  119  beam  paths 
and  two  principal  components  for  the  three  transitions  under  investigation.  Again,  the 
different  temperature  sensitivities  are  readily  apparent  in  the  absorption  maps.  Due  to  the 
low  ground  state  energy  of  the  225  cm'*  transition,  the  population  of  the  initial  state 
decreases  in  the  flame  due  to  thermal  excitation.  This  causes  the  absorption  to  decrease 
in  the  flame  resulting  in  an  absorption  map  with  a  “bowl”  in  the  center.  The  thermal 
excitation  from  the  flame  increases  the  initial  state  populations  of  the  1045  cm  '  and  1216 
cm  '  transitions  resulting  in  top-hat  shaped  absorption  maps  shown  here. 
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225  cm'^  Transition 


-3 


Figure  27. a.  Fit  absorption  map  for  the  225  cm”'  transition 


.  1045  cm’^  Transtition 

X  10 


Figure  27.b.  Fit  absorption  map  for  the  1045  cm''  transition 
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1216  cm‘^  Transition 


-3 


X  10 


0  0 

y  position  (cm)  ^ 

Figure  21. c.  Fit  absorption  map  for  the  1216  cm  *  transition 

Figure  27.  Example  frequency  integrated  absorption  coefficient  maps  for  each  of  the 

three  transitions  of  interest 

Combining  the  information  from  the  three  absorption  maps  using  a  Boltzmann 
plot  at  each  pixel  results  in  the  temperature  map  shown  in  Figure  28.  Again,  solving 
Equation  (22)  for  Nj  results  in  the  water  density  map  shown  in  Figure  28.  Using  the 
scaled  CFD  images  as  a  basis  for  the  lines  of  sight  allows  accurate  profile  determination 
for  the  flat  flame  burner. 
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17  Density  Map 
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Figure  28.  Composite  temperature  and  density  contour  plots  for  the  flat  flame  burner 
using  the  first  two  principal  components  and  119  beam  paths 

Taking  a  slice  from  the  two-zone  and  CFD  temperature  and  density  plots  allows  a 

comparison  to  the  thermocouple  data  of  the  flat  flame  burner.  Figure  29  shows  this 

comparison.  Constraining  the  cold  region  in  the  two-zone  model  to  agree  with  one  of  the 
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lines  of  sight  that  did  not  pass  through  the  flame  resulted  in  a  better  temperature  and 
density  determination  in  the  hot  region.  Both  the  CFD  and  the  two-zone  model,  however, 
give  very  good  agreement  on  mole  fraction  in  the  flame.  The  fit  compensated  for  the 
high  temperature  determination  of  the  CFD  model  in  the  center  by  lowering  the  density. 
Once  more  transitions  are  added  to  the  automated  integration  routine,  the  Boltzmann 
plots  should  become  more  accurate.  This  should  result  in  better  temperature  agreement 
with  both  the  two-zone  model  and  the  thermocouple  measurements.  Also,  because  the 
mole  fraction  in  the  hot  region  agrees  with  the  mole  fraction  from  the  complete 
combustion  model  used  in  the  thermocouple  experiment,  calculated  densities  should 
show  agreement  as  well. 
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Distance  Across  Burner  at  6  mm  Height  (mm) 


Region 

PCA 

Two-Zone 

TC/Hum.  Meter 

Hot  (40mm 
avg) 

T=  1712  +  89  K, 

]V[.F.  =  0.15±.02 

T  =  1655  +  17  K, 

M.F.  =  0.1 55 +  .009 

T=  1638±10K 

M.F.  =  0.1547 

Cold 

T  =  302  +  19K, 

M.F.  =  0.0073  ±4e-4 

T  =  302  +  4  K, 

M.F.  =  0.0076  +  4e-4 

T  =300  K  (Thermocouple) 

M.F.  =  0.0061  +  4e-4  (probe 
@293.5±.2K) 

Figure  29.  Comparison  of  CFD  results  to  two-zone  results  and  thermocouple 

measurements 


The  strategy  developed  in  this  chapter  was  able  to  obtain  accurate  temperature 
and  density  measurements  for  the  flat  flame  burner.  The  next  step  is  to  see  how  valuable 
the  strategy  is  when  applied  to  a  more  complicated  flow.  In  Chapter  VI,  the  same 
strategy  and  experimental  setup  used  for  the  flat  flame  burner  is  employed  to  investigate 
the  flow  provided  by  the  supersonic  wind  tunnel  in  Research  Cell  18. 
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V  Isolator  Experiment 


Introduction 

This  chapter  serves  as  an  introduction  to  the  isolator  experiment  performed  in 
Research  Cell  18.  The  motivation  for  the  flat  flame  burner  experiment  was  to  provide 
confidence  in  the  ability  of  TDLAS  based  tomography.  Its  performance  lends  credibility 
to  the  strategy  developed  and  its  implementation  in  more  complicated  flows  like  that  of 
the  isolator.  Paralleling  the  flat  flame  burner  experiment,  a  priori  information  was 
provided  by  CFD  calculations  to  reduce  the  number  of  beam  paths  necessary  to  perform 
tomography.  However,  the  CFD  calculations  did  not  agree  with  the  TDLAS  data  on  the 
location  of  the  first  reflection  in  the  shock  train.  This  inconsistency  made  it  impossible 
to  perform  tomography  without  first  adjusting  the  CFD  calculations.  Thus,  an 
optimization,  like  that  of  the  flat  flame  burner  diameter  and  center  location,  was 
performed  on  the  isolator  CFD  calculations  to  provide  a  better  match  to  the  TDLAS  data. 

With  16  lines  of  sight,  TDLAS  was  able  to  obtain  important  information  about  in¬ 
flow  static  pressure,  temperature,  and  density  before  any  tomography  was  performed. 
Vector  quantization  allowed  similar  conditions  to  be  averaged  to  provide  insight  into  the 
shock  structure.  This  information  enabled  comparison  with  the  CFD  and  a  similar 
parameterization  and  optimization  of  the  CFD  basis  images  paralleling  the  two-zone 
techniques  used  to  find  diameter  and  center  position  of  the  flat  flame  burner.  [16] 

Isolator  Experimental  Setup 

An  identical  data  collection  and  laser  control  system  to  the  system  shown  in 
Figure  11  was  used  to  make  measurements  in  Research  Cell  18.  However,  the 
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collimators  and  detectors  for  the  laser  beams  were  configured  in  metal  boxes  purged  with 
nitrogen.  This  eliminated  the  water  in  the  path  up  to  the  isolator  window  and  allowed 
measurements  to  be  made  only  in  the  region  of  interest.  The  isolator  provided  optical 
access  to  the  vitiated,  supersonic  flow  through  small,  approximately  1.5  by  3.5  inch, 
windows  on  each  side.  Methane  was  used  to  vitiate  the  flow,  providing  a  controlled 
amount  of  water  to  observe  using  TDLAS.  The  backpressure  of  the  flow  was  controlled 
by  a  valve  downstream  of  the  isolator.  Adjusting  the  backpressure  varied  the  position  of 
the  shock  in  the  isolator  which  was  then  observed  by  the  laser  array.  A  diagram  of  the 
wind  tunnel  is  shown  in  Figure  30.  Note  that  the  combustor  in  this  image  was  replaced 
with  a  ball  valve  to  provide  the  backpressure  resulting  in  a  controlled  shock  train. 


Figure  30.  Schematic  of  the  wind  tunnel  in  RC  1 8  at  Wright  Patterson  Air  Force  Base, 

adapted  from  reference[23] 

Closing  the  valve  increases  the  backpressure  which  pushes  the  shock  train  into  and 
through  the  isolator  where  the  laser  array  was  able  to  make  absorption  measurements. 
The  following  plots  show  the  CFD  model  of  the  entire  isolator  as  the  backpressure  is 
increased. 
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Figure  3 1 .  CFD  model  of  shock  in  the  isolator  as  backpressure  is  increased[3] 

The  image  in  Figure  32  shows  the  setup  of  the  diode  laser  array.  Note  that  this  is  a  very 
similar  setup  to  that  used  in  the  flat  flame  burner  experiment,  except  the  burner 
experiment  used  only  one  laser  array  and  one  photo-diode  detector  array.  For  the  isolator 
experiment,  eight  beams  were  used  on  each  side  giving  a  total  of  16  lines  of  sight.  Two 
of  the  beams  were  aimed  diagonally  across  the  flow  to  give  Doppler  information  and,  like 
the  angled  beams  in  the  flat  flame  experiment,  to  give  tomographic  information  in  the 
direction  perpendicular  to  the  flow. 
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Figure  32.  Image  of  experimental  setup  of  laser  diode  and  photo  diode  detector  arrays  for 

the  isolator  experiment[16] 


Comparison  to  CFD 

By  integrating  every  spectrum  at  the  data  collection  rate  of  1000  Hz,  instead  of 
averaging  over  the  entire  collection  like  in  the  flat  flame  burner  experiment,  oscillations 
in  the  shock  were  observed.  The  left  side  of  Figure  33  shows  a  waterfall  plot  of  the 
absorption  for  the  14  straight  lines  of  sight  as  time  progressed  at  a  constant  backpressure 
provided  by  the  ball  valve.  The  vector  quantization  algorithm,  originally  intended  to  aid 
in  training  set  dimensionality  reduction,  was  used  here  to  find  similarly  shaped  plots  of 
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the  absorption  across  the  window.  By  taking  the  average  of  these  similar  plots,  the  shock 
front  and  first  reflection,  which  are  the  first  and  second  pressure  peaks  in  the  shock  train, 
were  resolved  as  shown  in  the  two  plots  on  the  right  of  Figure  33. [16] 
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Figure  33.  Waterfall  plot  and  vector  quantized  averages  showing  shock  structure  in  the 

isolator 


Converting  the  plots  in  Figure  33  to  pressure  and  temperature  allows  a  comparison  with 
the  CFD  model.  The  plots  in  Figure  34  show  that  the  CFD  is  missing  the  first  reflection 
in  the  shock  train.  Using  the  data  from  only  14  lines  of  sight  perpendicular  to  the  flow, 
this  discrepancy  in  the  CFD  is  readily  apparent.  Even  without  tomography,  enough 
contour  information  is  available  from  the  TDLAS  to  make  a  comparison  with  the  CFD 
calculations. 

Performing  tomography  with  a  limited  number  of  beam  paths  relies  on  contour 
information  obtained  from  the  CFD  calculations.  The  discrepancy  between  the  CFD  and 
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line  of  sight  temperatures  and  pressures  shown  in  Figure  34  shows  that  the  current  CFD 
calculations  will  not  provide  an  adequate  basis  function  to  fit  the  data.  Tomography  was, 
therefore,  not  possible  using  the  given  CFD. 


Figure  34.  Averaged  absorption  plots  converted  to  temperature  and  pressure 


Optimization  of  CFD 

By  parameterizing  the  flow  in  terms  of  wall  thickness  and  shock  location,  an 
optimization  of  the  CFD  was  performed  to  correct  for  the  missing  first  reflection.  Similar 
to  the  scaling  of  the  diameter  of  the  CFD  for  the  flat  flame  burner,  the  position  of  the  first 
reflection  was  scaled  in  the  isolator  CFD  to  provide  a  better  fit  to  the  data  as  shown  in 
Figure  35. [16] 
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Figure  35.  Results  of  parameterization  and  optimization  of  CFD  images  to  TDLAS  data 
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The  top  two  images  show  the  top  down  view  of  the  original  CFD  temperature  and 
pressure  maps  for  the  isolator.  These  are  the  two-dimensional  representations  of  the  plots 
in  Figure  34.  Again,  the  first  reflection  is  missing  in  both  the  temperature  and  pressure 
maps.  After  the  optimization,  the  lowest  residual  was  achieved  by  reducing  thickness  of 
the  shock  and  moving  the  position  of  the  first  reflection  up  into  the  window.  Enough 
information  was  available  from  the  16  lines  of  sight  to  enable  a  modification  of  the  CFD 
to  better  fit  the  data.  Before  this  experiment,  CFD  calculations  of  this  sort  did  not  have 
any  in-flow  data  with  which  to  compare.  With  feedback  from  both  sides,  new  CFD 
calculations  can  be  made,  enabling  better  tomography  which  can  again  refine  the  CFD. 
Further  implementation  of  the  strategy  presented  in  this  thesis  will  put  the 
characterization  of  complicated  flows  with  limited  optical  access  in  reach. 
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VI  Conclusions  and  Recommendations 


The  primary  objective  of  this  thesis  was  to  provide  validation  of  the  developed 
strategy  to  perform  tomography  using  a  limited  number  of  TDLAS  beam  paths.  By 
thoroughly  characterizing  the  flat  flame  burner  through  the  traditional  methods  of 
thermocouples,  pressure  probes,  and  humidity  probes,  a  good  point  of  comparison  was 
obtained.  Using  fairly  basic  geometry,  the  two-zone  analytic  method  for  determining 
temperature  of  the  hot  and  cool  region  showed  very  good  agreement  with  the 
thermocouple  data.  Applying  this  same  model  to  perform  tomography  with  varying 
numbers  of  beam  paths  demonstrated  not  only  good  agreement  with  both  the  two-zone 
analytic  method  and  thermocouple  method,  but  that  good  agreements  can  be  achieved 
with  very  limited  optical  access.  Finally,  performing  tomography  with  the  CFD  model 
allowed  excellent  determination  of  the  contour  of  the  burner  profde  as  well  as  the  mole 
fraction  of  water  in  the  flame.  By  stepping  through  each  level  of  tomography,  this  thesis 
shows  that  tomography  of  a  flow  with  limited  optical  access  is  possible  using  TDLAS.  It 
provides  a  proof  of  concept  for  the  more  geometrically  challenging  isolator  experiment. 

The  utility  of  both  the  VQ  and  PCA  algorithms  was  tested  for  tomography.  VQ 
proved  to  be  helpful  in  the  near  real-time  application  of  tomography  in  the  isolator.  It 
worked  very  well  as  a  method  to  optimize  parameters  in  both  the  two-zone  model  and  the 
CFD  model.  It  was  also  useful  in  selecting  similarly  shaped  shock  structures  in  the  time 
resolved  shock  oscillations  of  the  isolator  flow.  PCA  was  able  to  reduce  the  CFD 
training  sets  used  for  the  flat  flame  burner  from  175  to  just  3  images,  which  drastically 
reduced  the  number  of  beam  paths  necessary  to  perform  tomography. 
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This  thesis  also  demonstrated  the  technology  necessary  to  perform  the  TDLAS 
measurements.  Using  a  relatively  inexpensive  system,  up  to  16  beam  paths  were  sampled 
over  three  spectral  regions  at  sweep  rates  greater  than  1  kHz  resulting  in  48  spectra  per 
sweep.  This  high  data  rate  allowed  near  real-time  reconstruction  of  the  flow  in  both  the 
flat  flame  burner  as  well  as  the  isolator. 

The  flat  flame  burner  results  lay  the  groundwork  for  understanding  the  effects  of 
temperature  changes  on  TDLAS  based  tomography.  The  method  developed  for  obtaining 
tomographic  maps  of  temperature  and  density  has  been  carried  over  to  the  isolator 
experiment,  where  the  effects  of  shock  train  structure  on  TDLAS  based  mass- flux  sensors 
were  studied.  The  results  of  the  isolator  experiment,  introduced  in  Chapter  5,  show  that 
this  approach  can  be  carried  over  to  other  parts  of  the  supersonic  combustion  flow  path 
that  involve  complex  flow  distributions,  such  as  the  exhaust  plane  of  supersonic 
combustors.  [16] 

Recommendations  for  future  work  are  manifold.  In  the  short  term,  a  more 
rigorous  method  of  automating  the  integration  of  absorption  spectra  is  necessary. 
Providing  more  points  for  the  Boltzmann  plots  should  increase  the  accuracy  of  the 
temperature  determination  because  now,  only  three  points  are  used  to  perform  the  linear 
fit.  Error  in  just  one  point  can  cause  a  significant  difference  in  the  temperature 
determination.  With  more  points,  the  error  will  be  reduced.  Also,  the  current  method  for 
determining  absorption  coefficient  maps  fits  each  transition  independently.  Performing 
the  fit  with  all  of  the  transitions  at  once  will  help  to  constrain  the  fit  and  provide  a  more 
accurate  result.  A  new  flat  flame  burner  with  four  independent  quadrants  is  currently 
being  set  up.  Performing  the  same  analysis  on  this  more  complicated  flow  would  help  to 
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perfect  the  method  and  provide  further  proof  of  concept.  Future  work  could  also  include 
helping  to  refine  the  CFD  model  by  iterating  on  the  testing  of  the  CFD  and  repeating  it  to 
get  a  better  fit. 

Scheduled  future  work  also  includes  flight  testing  of  a  TDLAS  system  on  board  a 
SCRAM  Jet  engine.  This  will  provide  a  great  deal  of  insight  into  how  well  TDLAS  can 
perform  in  the  harsh  environments  encountered  by  a  supersonic  vehicle.  TDLAS  has  the 
potential  to  revolutionize  the  field  of  engine  diagnostics,  particularly  in  flight  vehicles 
where  weight  and  size  are  of  critical  importance. 
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Appendix  A.  Beer’s  Law  to  Frequency  Integrated  Absorption  Coefficient 


Adapted  from  Bemath  pp  272-274. [8] 
Beer’s  Law: 


I  =  I^e  =I,e 

Assuming  thermodynamic  equilibrium,  substitute  in  the  following  equations: 

_  §l  ^-hviJkT 

No  go 
-'''o  ~ 

q 


Then, 


,-^0  IkT 


So  gl  -hv]0  I  kT 


I,e 


Si  So 


)/ 


I,e 


^SO^T^  ^  ^-hvipIkT  ^ 


Rearrange: 


■In 


\^oJ 


=  a 


g^Nje 


-E^/kT 


(1- 


-hv^Q/kT 


Substitute  in  for  the  cross  section: 


(7  = 


3e^hc  go 


And  get: 


In 


\^oj 


3^o^c  go  q 


-E^/kT 


.(l_e-»''^o/kTy 


But  this  can  be  simplified  using  the  HITRAN  line  intensity: 

2;rV,n5'„ 


S'=- 


10^10  ^-E^/kT 


Se^hcq 


So, 
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■In 


=  S\T)-Nj-l-g{v-v,^)=A  =  a-l 


Integrate  over  all  frequency  to  eliminate  the  line  function  and  get: 

p  00  p  00 

[  A{v)dv=\  -In 

J  -00  J  -00 

Where  a  ■  I  has  units  of  cm  '  if  v  is  in  wavenumbers. 


^  l"\ 

—  dv  =  S'NJ  =  a-l 

J  ‘ 

\^oJ 
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Appendix  B.  CFD  to  Frequency  Integrated  Absorption  Coefficient  Map 

For  the  flat  flame  burner,  the  CFD  calculations  give  temperature  in  K,  mass  fraction  of 

H2O,  and  the  total  density  of  the  burning  gas  in  kg/m  all  as  a  function  of  radial  position. 

First,  the  total  gas  density,  D^?,  and  the  mass  fraction  of  water,  Fhio,  are  used  to  calculate 

the  density  of  water  in  the  flame  in  molecules/cm^,  Nf. 

p,  1  amu  IH2O  molecule  Im^ 

^“'1.66x10-"’ kg  18.0152  amu  lO*"  cm' 


•  3.34  X  10» 


Next, 


a  =  S\T)N^  = 


3hcq{T^^f) 


/kT^ef 


_  ^-hcVi„lkT  ^ 


q{T) 


Eq  I 


(1-e- 


■N^ 


ref  ' 


=  S'iT„,) 


•ef )  ^ 


-E^IkT 


q{T)  _ ^“"^'''10 ' '‘V./  ^ 


-Acvi  n  /  kT^^f  N  T 


where  g'C^re/)  has  units  of  cm/molecule  and  Nt  has  units  of  molecules/cm  ,  giving  a 

units  of  cm'^.  The  partition  functions,  q(Tref)  and  q(T)  are  calculated  using  the  Matlab 
function  called  partition  (found  in  Appendix  C)  written  by  Dr.  Chad  Lindstrom.  The  line 
intensity,  *S"(^re/) » units  of  cm/molecule,  the  term  energy,  Eq,  and  the  transition  energy. 


V]o,  in  units  of  wavenumbers,  are  given  by  a  HITRAN  table. 
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Appendix  C.  Matlab  Programs 


The  following  programs  were  written  by  Dr.  Chad  Lindstrom  and  Ryan  Givens  to 
ereate  basis  funetions  for  the  flat  flame  burner  experiment,  fit  the  basis  functions  to  the 
data,  extract  temperature  and  density,  and  analyze  the  results: 

Burnerbasisv2: 

function  [A^B^x^y]  = 

burnerbasis (C_list , density, temp, r ratio, params , prof ilex, maxT) 

%  [B,x,y]  = 

%  burnerbasisv2 (C_list , density, temp, r ratio, params , prof ilex, maxT) 

%  function  to  create  snapshots  for  a  McKenna  burner 

o 

o 

%  OUTPUT: 

%  B  =  a  PxQxM  tensor  of  the  M  = 

%  numelements (C_list) ^numelements (D_list) snapshots 
%  X  =  The  X  coordinates  of  the  matrix 
%  y  =  The  y  coordinates  of  the  matrix 

o  o 
o  o 

%  INPUT: 

%  C_list  =  A  Nx2  vector  (xO  yO;xf  yf;dx  dy)  that  describes  the  range 
%  that  the  center  of  the  burner  will  be  moved.  It  is  on  a  square  grid 
%  that  goes  from  xO  to  xf  and  yO  to  yf  in  increments  of  dx  and  dy 

o 

o 

%  density  =  A  NxM  matrix  consisting  of  a  the  mass  fraction  and  density 
%  of  the  various  basis  guesses  in  the  form  (MFl  D1  MF2  D2  . . . ) 

o 

o 

%  temp  =  A  Nx(M/2)  matrix  of  the  radial  temperature  profiles 
%  rratio  =  A  list  of  ratios  for  the  max  water  density  in  the  flame  to 
%  the  min  water  density  in  the  flame 

%  params  =  A  7x1  vector  with  (Lx, Ly, dx, dy, Egnd, Sref , nu)  where  Lx  is  the 
%  length  of  the  region,  Ly  is  the  length  of  the  region  in  the  y- 
%  direction,  dx  is  the  distance  between  adjacent  pixels,  dy  is  the 
%  distance  between  the  adjacent  pixels  in  the  y-direction,  Egnd  is  the 
%  ground  state  energy,  Sref  is  the  line  intensity  at  296  K,  and  nu  is 
%  the  wavenumber  of  the  transition 

%  maxT  =  A  list  of  maximum  temperatures  (RT  is  set  to  293  K  =  20  deg  C) 

o 

o 

%  Author:  Chad  Lindstrom 

Nx  =  floor ( (C_list (2, 1) -C_list (1, 1) ) /C_list (3, 1) ) 

Ny  =  floor ( (C_list (2, 2) -C_list (1, 2) ) /C_list (3, 2) ) 

NB  =  size (temp, 2 ) ;  %  number  of  basis  functions 
No  =  size  (rratio, 1 ) ; 

Nt  =  size (maxT, 1 ) ; 

index2  =  1; 

for  i=l:NB  %calculate  ID  absorption  profiles  for  all  cases 
for  j  =1 : No 

for  k=l : Nt 

D  =  density  (:,  2 "^i-l density  (:,  2 "^i )  "^3 . 34 4el 9 ; 
maxD  =  D ( 1 ) ; 
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minD  =  maxD/rratio ( j ) ; 

D  =  scaleshif t (maxD, minD, D) ; 

T  =  scaleshif t (maxT ( k) , 2 92 , temp ( : , i ) ) ; 

A(:,index2)  =  alpha (T, D, params (5) , params ( 6) , params (7 ) ) 

index2  =  index2+l; 

end 

end 

end 

index2=index2-l  ; 

center  =  ones (2,1); 
center (1)  =  C_list(l,l); 
center (2)  =  C_list(l,2); 

[H,x,y]  =  burnerv2 ( center , A (:, 1 ), prof ilex, params ) ; 

LI  =  size (H, 1) ; 

L2  =  size (H, 2) ; 

B  =  ones  (Li, L2,Nx^Ny^NB^No)  ; 
index  =  1 ; 

for  i  =  1 : Nx 

for  j  =1 : Ny 

for  k=l : index2 

center  (1)  =  C_list  ( 1 , 1 )  +  ( i-1 ) '^C_list  ( 3 , 1 )  ; 
center(2)  =  C_list ( 1 , 2 ) + ( j -1 ) ^C_list ( 3 , 2 ) ; 

[H,x,y]  =  burnerv2 ( center , A (:, k) , prof ilex, params ) ; 
B ( : , : , index)  =  H; 
index  =  index+1; 
end 

end 

end 

end 
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Alpha: 

function  A  =  alpha (T, N, Egnd, Sref, nu) 

%  function  A  =  alpha (T, N, Egnd, Sref, nu) 

%  is  a  function  that  calculates  the  integrated  absorption  coefficient 
%  map  given  the  temperature  and  density  maps  as  well  as  the  parameters 
%  for  the  specific  spectroscopic  transisition 

o 

o 

%  Outputs: 

o 

o 

%  A  =  the  absorption  map  in  cm^-1 

o 

o 

%  Inputs: 

o 

o 

%  T  =  The  temperature  in  Kelvin  (may  be  a  vector  or  matrix) 

%  the  maximum  temperature  supported  is  2990  K  and  the  minimum 

%  is  200  K 

%  N  =  The  number  density  of  the  molecules  (cm'^-S)  must  be  the 

%  same  size  as  T 

%  Egnd  =  The  ground  state  energy  in  cm^-1  that  is  being  probed 

%  Sref  =  The  line  intensity  from  Hitran  at  296  K 

%  nu  =  wavenumber  of  exciting  light  (cm"^-l) 

%  Author:  Chad  Lindstrom 

Qref  =  partition (2 96, nu) ; 

Q  =  partition (T, nu) ; 

S  =  Sref^Qref^exp(-1.439839851^Egnd^ (1./T-1/300)  )  ./Q; 

A  =  N.^S; 
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Partition: 


function  Q  =  Partition (T, nu) 

o 

o 

%  function  Q  =  Partition (T, nu) 

%  is  a  function  that  calculates  the  partition  function  and  accounts  for 
%  population  in  the  excited  state  based  on  the  Hitran  parameters 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 


Outputs : 

Q  = 


Inputs : 


The  value  of  the  partition  function  at  temperature  T  (if  T 
is  a  vector  or  matrix  Q  will  be  the  same  size  vector  or 
matrix) 


%  T  =  The  temperature  in  Kelvin  (may  be  a  vector  or  matrix) 

%  the  maximum  temperature  supported  is  2990  K  and  the  minimum 

%  is  200  K 

o 

o 

%  nu  =  wavenumber  of  exciting  light  (cm"^-l) 

o 

o 

%  Author:  Chad  Lindstrom 


R1  =  [-1.822185^10  4.230752^10^-1  7.703331^10^-4  6.185911^10^-9]; 

R2  =  [-6.594415^10  6.945487^10^-1  2.468324^10^-4  3.484234^10^-7]; 

R3  =  [-6.18536^10  6.973356^10^-1  2.183351^10^-4  3.699410^10^-7]; 

R4  =  [1.993301^10^5  -4.280657^10^2  3.073247^10^-1  -7.289637^10^-5]; 
R5  =  [-1.311914^10^3  3.214013  -1.522213^10^-3  7.885532^10^-7]; 

R6  =  [-4.529464^10'^3  7.314166  -3 . 2  6518  9^10'^-3  1 . 035648^10'^-6]  ; 


rows  =  size (T, 1 ) ; 
columns  =  size  (T, 2); 


for  i  =  1 : rows 

for  j  =  1: columns 
if  (T(i, j)<=500) 

Q(i,j)  =  Rl(l)+Rl(2)^T(i,  j)+Rl(3)^T(i,  j)"2+Rl(4)^T(i,  j)"3; 
elseif  (T  (i, j ) <1000) 

Q(i,j)  =  R2(l)+R2(2)^T(i,j)+R2(3)^T(i,j)-2+R2(4)^T(i,j)-3; 
elseif  (T  (i, j ) <1350) 

Q(i,j)  =  R3(l)+R3(2)^T(i,  j)+R3(3)^T(i,  j)'^2+R3(4)^T(i,  j)'^3; 
elseif  (T  (i, j ) <1425) 

Q (i, j )  =  R4 (1) +R4 (2) ^T (i, j ) +R4 (3) ^T (i, j ) ^2+R4 ( 4 ) ^T ( i , j ) ^3 ; 
elseif  (T  (i, j ) <2500) 

Q(i,j)  =  R5(l)  +  R5(2)^T(i,  j)+R5(3)^T(i,  j)'^2+R5(4)^T(i,  j)'^3; 
else  %if  (T  (i, j)<2990) 

Q(i,j)  =  R6(l)  +  R6(2)^T(i,  j)+R6(3)^T(i,  j)"2+R6(4)^T(i,  j)"3; 

end 

Q(i,j)  =  Q(i,j)/(l-exp(-4.1357^3^nu/(8.617^T(i,j)))); 
end 

end 
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Scaleshift: 

function  Pn  =  scaleshift ( PO , Pmx, data) 

%  function  Pn  =  scaleshift ( PO , Pmx, data) 

%  is  a  function  that  will  scale  and  shift  an  array  for  the  burner 
%  problem  so  that  the  temperatures  or  concentrations  match  measured 
%  values  at  the  center  and  edges  of  the  profiles 

o 

o 

%  Outputs: 

o 

o 

%  Pn  =  the  new  profile 

o 

o 

%  Inputs: 

o 

o 

%  PO  =  The  profile  value  at  the  center  of  the  burner 

%  Pmx  =  The  profile  value  at  the  edge  of  the  burner  region 

%  data  =  The  profile (s)  (each  column  is  taken  to  be  a  profile) 

o 

o 

%  Author:  Chad  Lindstrom 

columns  =  size (data, 2 ) ; 
length  =  size (data, 1 ) ; 

for  i=l: columns 

xO  =  data ( 1 , i) ; 
xf  =  data ( length, i ) ; 
m  =  (Pmx-PO) / (xf-xO)  ; 
b  =  PO-m'^xO; 

Pn  =  data ( : , i) ^m+b; 

end 
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bumerv2: 

function  [burnerprof ile, x, y]  =  f (center, prof ile, prof ilex, params) 

%  [burnerprof ile,  X,  y]  =  f (center, prof ile, prof ilex, params)  function  to 
%  create  snapshots  for  a  McKenna  burner 

o 

o 

%  OUTPUT: 

%  burnerprof ile  =  a  profile  that  has  an  inner  circle  with  magnitude 
%  ratio  and  an  outer  rectangular  region  with  a  value  of  1 

o 

o 

%  K  =  The  X  coordinates  of  the  matrix 

%  Y  =  The  y  coordinates  of  the  matrix 

o  o 
o  o 

%  INPUT: 

%  center  =  A  2x1  vector  (xO  yO)  that  describes  the  center  of  the  burner 
%  profile  =  The  ID  radial  profile  for  the  burner  (assumes  last  point  is 
%  absorbence  coefficient  in  ambient  air) 

%  profilex  =  The  x  coordinates  for  the  burner  profile 

%  params  =  A  4x1  vector  with  (Lx, Ly, dx, dy)  where  Lx  is  the  length  of 

%  the  region,  Ly  is  the  length  of  the  region  in  the  y-direction,  dx  is 

%  the  distance  between  adjacent  pixels,  dy  is  the  distance  between  the 
%  adjacent  pixels  in  the  y-direction 

o 

o 

%  Author:  Chad  Lindstrom 


Nx  =  ceil (params ( 1 ) /params ( 3 )) +1 ; 

Ny  =  ceil (params (2 ) /params ( 4 )) +1 ; 

X  =  ( 0  :  Nx-1 )  "^params  ( 3 )  ; 
y  =  ( 0  :  Ny-1 )  "^params  ( 4  )  ; 

length  =  size  (prof ile, 1 ) ; 

burnerprof  ile  =  ones  (Nx,  Ny)  "^prof  ile  ( length)  ; 
centerpixelx  =  center ( 1 ) /params ( 3 ) ; 
centerpixely  =  center (2) /params (4) ; 

Rpixelx  =  prof ilex ( length) /params ( 3 ) ; 

Rpixely  =  prof ilex ( length) /params ( 4 ) ; 
minpixelx  =  floor ( centerpixelx-Rpixelx) ; 
maxpixelx  =  ceil (centerpixelx+Rpixelx) ; 
minpixely  =  floor ( centerpixely-Rpixely) ; 
maxpixely  =  ceil (centerpixely+Rpixely) ; 

% - Fill  the  snapshot  matrix - 

for  i  =  minpixelx : maxpixelx 

for  j  =  minpixely : maxpixely 

h  =  sqrt((x(i)-center(l))^2+(y(j)-center(2))^2); 

%burnerprof ile  ( i ,  j  )  =  erf c  (  (h-D/2 )  )  "^  (ratio-1 ) /2  +  1  ; 

if  (h<prof ilex ( length) ) 
if  (h>prof ilex ( 1 ) ) 

burnerprof ile ( i ,  j )  =  interplq (prof ilex, prof ile, h) ; 
else 

burnerprof ile ( i , j )  =  interplq (prof ilex, prof ile, prof ilex (2 )) ; 
end 

end 

end 

end 

burnerprof ile  =  burnerprof ile ' ; 
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Create_basisv2: 

function  [basis , eig_val ]  =  f (snapshots, mval, N) 

%  [basis , eig_val ]  =  Create_basis (snapshots, mval, N)  function  to  create 
%  basis  from  snapshots  for  large  data  sets  using  SVD 

o 

o 

%  OUTPUT: 

%  basis  =  a  PxQxM  tensor  of  M  basis  functions 
%  eig_val  =  an  increasing  vector  of  the  eigenvalues 

o 

o 

%  INPUT: 

%  snapshots  =  a  PxQxM  tensor  of  the  M  snapshots 
%  mval  =  the  mean  value  used  to  subract  from  snapshots 
%  N  =  number  of  basis  elements  to  retain 

o 

o 

%  Author:  Chad  Lindstrom 

M  =  size  ( snapshots , 3 ) ;  %Determine  number  of  snapshots 
rows  =  size  ( snapshots , 1 ) ; 

mvalv  =  matrixtovector (mval ) ; 
sv  =  matrixtovector ( snapshots ) ; 
for  i  =  1 : M 

c(i)  =  mvalvXsv ( :  ,  i)  ; 

sv(:,i)  =  sv ( : , i ) -mvalv^c ( i ) ;  %subract  mean 

end 

clear  mvalv; 

% - Compute  elements  of  Matrix  to  Determine  Eigenvalues - 

sv  =  sv/ sqrt (M) ; 

[U,S,V]  =  svd(sv,0); 
clear  sv; 

eig_val  =  diag(S.^2); 
eig_val  =  f lipud (eig_val) ; 
for  i  =  1 : N 

basis (:,:,!)  =  snapshots (:,:,!); 
basis (:,:,!)  =  0  ; 
for  j  =  1 : M 

basis  (:,:,!)  =  V  ( j  ,  i)  (snapshots  (  :  ,  :  ^  j  )  - 
c  ( j  )  "^mval)  +basis  (:,:,!); 
end 

end 

basis  =  f lipdim (basis , 3 ) ; 
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Matrixtovector : 

function  matrixtovector  =  f (matrix) 

%  matrixtovector  =  f (matrix)  function  to  create  vector  from  matrix 

o 

o 

o 

o 

%  OUTPUT: 

%  matrixtovector  =  either  a  single  vector  comprised  of  each  column  of 
%  the  matrix  stacked  on  top  of  each  other  if  only  one  input  matrix  is 
%  given  or  a  matrix  with  each  column  vector  corresponding  to  each 
%  matrix  of  the  3rd  rank  tensor 

o 

o 

%  INPUT: 

%  matrix  =  either  a  PxQ  matrix  or  a  PxQxM  tensor  with  M  PxQ  matrices 

o 

o 

%  Author:  Chad  Lindstrom 


M  =  size (matrix, 3 ) ;  %Determine  size  of  matrix 
P  =  size (matrix, 1 ) ; 

Q  =  size (matrix, 2 ) ; 
matrixtovector  =  ones  ( P'^Q, M)  ; 

% - Change  matrix  (es)  to  vector  (s) - 

for  i  =  1:M 
for  j  =1:Q 

row_start  =  (j-l)'^P+l; 
row_end  =  j'^P; 

matrixtovector (row_start : row_end, i) ^matrix ( : , j , i) ; 

end 

end 
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Datafit2: 

Function  [bestfit^resd, fiterr] =dataf it (mval, basis, data, x, y, Beams, sigma) 
%  [bestf it , resd, fiterr ]  =calculatef it (mval, basis, image, x, y, Beams, sigma) 
%  Calculates  the  fit  to  a  real  data  vector 

o 

o 

%  OUTPUT: 

o 

o 

%  bestfit  =  the  bestfit  for  the  data 
%  resd  =  the  fit  residuals 

o 

o 

%  Input: 

o 

o 

%  mval  =  The  mean  matrix 
%  basis  =  The  basis  used  for  the  fit 
%  data  =  The  data  vector 
%  X  =  X  coordinates  of  the  image 

%  y  =  y  coordinates  of  the  image 

%  Beams  =  The  beams  matrix  that  describe  the  beams  for  inversion 

%  sigma  =  The  error  in  the  estimation  of  the  fit  integrals 

o 

o 

%  Author:  Chad  Lindstrom 

[condnum, mval_pro j ]  =  Cond_Proj(Beams,mval,x,y,300); 

[condnum, basis_pro j ]  =  Cond_Proj  (Beams , basis , x, y, 300 ) ; 

%V  =  diag  (ones  ( 1 ,  size  (data,  1 ))  ^sigma'^P )  ; 

V  =  diag ( sigma . ^2 ) ; 

[cm, stdcm, mse]  =  Iscov (mval_pro j , data, V) ; 

stdcm  =  stdcm^sqrt ( 1/mse) ; 

diff_proj  =  data-cm^mval_pro j ; 

bestfit  =  cm'^mval; 

fiterr  =  ( stdcm'^mval_pro j  )  .^2; 

V  =  diag ( fiterr ) +V; 

[c,stdc,mse]  =lscov (basis_pro j , dif f_pro j , V) ; 
stdc  =  stdc^sqrt ( 1/mse) ; 

N  =  size  (c) ; 
bestfit  =  cm'^mval; 
fiterr  =  ( stdcm^mval ) . ^2 ; 
for  i=l:N 

bestfit  =  bestfit  +  c  (i)  "^basis  (  :  ,  :  ,  i)  ; 
fiterr  =  f  iterr+ ( stdc  ( i )  "^basis  (  :  ,  :  ,  i )  )  .  ^2  ; 

end 

resd  =  data  -  cm'^mval_pro  j -basis_pro  j  "^c; 
fiterr  =  sqrt ( fiterr ) ; 

return 
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Mapthattempne  w : 

function 

[tempmap, errmap, bolt zmannxs , bolt zmannys ] =mapthattempnew (al391_la, al391_ 

Ib,al391_2a,al391_2b,al391_2c,al391_2d,al397a,al397b,al397c,errl391_la7 

errl391_lb, errl391_2a, errl391_2b, errl391_2c, errl391_2d, errl397a,errl397 
b, errl397c) 

%  mapthattempnew  is  a  function  that  creates  a  temperature  map  and 
%  errormap  using  the  frequency  integrated  absorption  coefficients  and 
%  their  corresponding  errors  from  the  line  of  sight  data  using  a 
%  Boltzmann  Plot  at  each  pixel 

o 

o 

%  boltzmannxs  and  boltzmannys  are  the  x  and  y  data  points  for  the 
%  Boltzmann  plot 

o 

o 

%  Author:  Ryan  Givens 
k=.6950356; 

el391_la=1045.0579;  sl3 91_la=8 . 335E-22 ; 
el391_lb=1474 .9808;  sl391_lb=1.131E-23; 
el391_2a=136. 7617;  sl391_2a=l .488E-20; 
el391_2aa=1216. 1945;  sl391_2aa=7.273E-23; 
el391_2b=1216. 1945;  sl391_2b=2.438E-22; 
el391_2c=224 .8384;  sl391_2c=5.577E-22; 
el391_2d=1477 .2974;  sl391_2d=3.277-23; 
el397a=224 .8384;  sl397a=l .168e-20; 
el397b=222 .0528;  sl3 97b=4 . 2 12e-2 1 ; 
el397c=1874 .973;  sl397c=5 . 983e-24; 

x_points= [el391_la/k,  ... 

. . .el391_lb/k,  . . . 

.  . .el391_2a/k,  . . . 

el391_2b/k,  ... 

.  .  .  el391_2c/k,  .  .  . 

.  .  .el391_2d/k,  .  .  . 

el397a/k,  .  .  . 

...el397b/k,  ... 

. . .el397c/k  . . . 

]  ; 


vl391_la=7185. 59731; 
vl391_lb=7185. 40077; 
vl391_2a=7181 . 15578; 
vl391_2aa=7181 .09528; 
vl391_2b=7179.75201; 
vl391_2c=7180.39972; 
vl391_2d=7180 . 61301; 
vl397a=7161 .4101; 
vl397b=7160 .23543; 
vl397c=7160.8121; 

al391_2a=al391_2a- (sl391_2aa/sl391_2b) .^al391_2b;  %correct  for 
unresolved  line 


bolt z_mapl3 91_la=log  (al3  91_la/ sl3  91_la'^exp  (-el3  91_la/k/ 300)  (1-exp  (- 
vl391  la/k/300) ) ) ; 
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boltz_mapl391_lb=log (al391_lb/sl391_lb^exp (-el391_lb/k/300 ) ^ (1-exp ( 
vl391_lb/k/300) ) ) ; 

bolt z_mapl3 91_2a=log  (al391_2a/sl391_2a'^exp  (-el391_2a/k/300)  ^  (1-exp  ( 
vl391~2a/k/300) ) ) ; 

boltz_mapl391_2b=log (al391_2b/sl391_2b^exp (-el391_2b/k/300) ^ (1-exp ( 
vl391_2b/k/300) ) ) ; 

boltz_mapl3 91_2c=log  (al3  91_2c/ sl3  91_2c'^exp  (-el391_2c/k/ 300 )  ( 1-exp  ( 

vl391_2c/k/300) ) ) ; 

boltz_mapl391_2d=log (al391_2d/sl391_2d^exp (-el391_2d/k/300 ) ^ (1-exp ( 
vl391_2d/k/300) ) ) ; 

boltz_mapl3  97a=log  (al3  97a/ sl3  97a'^exp  (-el397a/k/300)  ( 1-exp  (- 

vl397a/k/300) ) ) ; 

boltz_mapl397b=log (al397b/sl397b^exp (-el397b/k/300) ^ (1-exp (- 
vl397b/k/300) ) ) ; 

boltz_mapl397c=log (al397c/sl397c^exp (-e 1397c/ k/ 300) ^ (1-exp (- 
vl397c/k/300) ) ) ; 

weight_mapl3 91_la=al3 91_la . /errl391_la; 
weight_mapl391_lb=al391_lb . /errl391_lb; 
weight_mapl3 91_2a=al3 91_2a . / errl391_2a; 
weight_mapl3 91_2b=al3 91_2b . /errl391_2b; 
weight_mapl3 91_2c=al3 91_2c . /errl391_2c; 
weight_mapl3  91_2d=al3  91_2d . /errl3  91_2d; 
weight_mapl397a=al397a ./errl397a; 
weight_mapl397b=al397b ./errl397b; 
weight_mapl3 97c=al3 97c ./errl397c; 

tempmap=ones (size(al397a) )  ; 
er rmap=tempmap ; 

bolt zmannxs=x_points ; 

boltzmannys=ones (size (al397a, 1) , size (al397a, 2) , length (x_points) ) ; 

for  i=l : size (al397a, 1 ) 

for  j=l : size (al397a,  2 ) 

y_points= [boltz_mapl3 91_la (i, j ) ,  ... 

. . . bolt z_mapl3 91_lb (i, j ) ,  ... 

.  .  . bolt z_mapl3 91_2a (i, j )  ,  ... 

boltz_mapl3  91_2b (i, j )  ,  ... 

. . . bolt z_mapl3 91_2c (i, j ) ,  ... 

. . .boltz_mapl391_2d (i, j ) ,  ... 

boltz_mapl397a (i,  j )  ,  ... 

.  .  . bolt z_mapl3 97b (i,  j  )  ,  ... 

.  . . boltz_mapl397c (i,  j  )  ... 

]  ; 

weights^ [weight_mapl3  91_la (i,  j  )  ,  ... 

. . . weigh t_mapl 39 l_lb (i, j ) ,  ... 

. . . weight_mapl3 91_2a (i, j ) ,  ... 

weight_mapl3 91_2b (i, j )  ,  ... 

. . . weight_mapl391_2c (i, j ) ,  ... 

. . . weight_mapl3 91_2d ( i , j ) ,  ... 

weight_mapl3  97a (i,  j )  ,  ... 

.  .  . weight_mapl3 97b (i, j )  ,  ... 

...  weigh t_mapl3 97c ( i ,  j  )  ,  ... 

]  ; 

S= [ones (length (x_points) , 1) ,  x_points ' ] ; 
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[p, se, sd] =lscov (S, y_points ' , weights ' ) 
tempmap (i, j ) =-1/ (p (2 ) ) ; 
errmap (i,j)=l/ (p(2) ) ^2^se(2)  ; 
boltzmannys (i, j  a  : ) =y_points; 

% tempmap (180,100) 

end 

end 
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Reperr2: 

function  [maxerr , avgerr , cmean]  =  reperr2 (B, basis , mval , nb) 

%  [maxerr , avgerr , cmean]  =  reperr2 (B, basis , mval , nb) 

o 

o 

o 

o 

%  OUTPUT: 

%  maxerr  =  the  maximum  absolute  percent  error  observed  for  all  elements 
%  of  the  training  set 

o 

o 

%  avgerr  =  the  average  absolute  percent  error  observed  for  all  elements 
%  of  the  training  set 

%  cmean  =  the  mean  value  of  the  expansion  coefficients 

o 

o 

%  INPUT: 

%  B  =  a  PxQxM  tensor  of  snapshots  (training  set) 

%  basis  =  a  PxQxM  tensor  of  basis  elements  created  from  B 
%  nb  =  the  number  of  basis  elements  to  retain  (the  largest 

o 

o 

%  Author:  Chad  Lindstrom 

Bv  =  matrixtovector (B) ; 
basisv  =  matrixtovector (basis ) ; 
mvalv  =  matrixtovector (mval ) ; 
lengthl  =  size  (basisv, 2 ) ; 
length2  =  size(Bv,2); 
cmean  =  ones ( lengthl , 1 ) ; 
cmean  =  0; 
avgerr  =  [0  0 ]  ; 
maxerr  =  [0  0 ]  ; 

for  i=l:length2 

cav  =  mvalv\Bv ( : , i ) ; 

Bv2  =  Bv(:,i)  -  cav^mvalv; 

c  =  basisv (:, lengthl-nb+1 : lengthl ) \Bv2 ; 
cmean  =  cmean+abs (c) ; 

result  =  basisv  (:,  lengthl-nb+1 :  lengthl )  "^c+cav^mvalv; 
abspererror  =  abs (Bv (:, 1) -result ) ^100 . /abs (Bv (:, 1) ) ; 
mr  =  max (abspererror )  ; 
if  (mr>maxerr ( 1 ) ) 

maxerr  =  [mr  1] ; 

end 

avr  =  mean (abspererror )  ; 
if  (avr>avgerr ( 1 ) ) 

avgerr  =  [avr  1] ; 

end 

end 

cmean  =  cmean/length2 ; 


95 


CreateVQbasis: 

function  [basis , which] =f ( snapshots , initial_codebook, N, epsilon) 

%  [basis , which]  =  Create_VQbasis ( snapshots , initial_codebook, N, epsilon) ; 
%  function  to  create  basis  from  snapshots  for  large  data 

o 

o 

%  OUTPUT: 

%  basis  =  a  PxQxM  tensor  of  M  basis  functions 

%  which  =  a  IXM  vector  indicating  which  codevector  each  training  set 
%  image  is  associated  with 

o 

o 

%  INPUT: 

%  snapshots  =  a  PxQxM  tensor  of  the  M  snapshots 
%  initial_codebook  =  PxQxN  tensor  of  initial  codevectors 
%  N  =  number  of  basis  elements  desired 

%  epsilon  =  distortion  measure  (small  number  by  which  to  adjust  average 
%  basis  vectors 

o 

o 

%  Author:  Ryan  Givens 


M  =  size  ( snapshots , 3 ) ;  %determines  the  number  of  snapshots  given 
T  =  size ( initial_codebook, 3 ) ;  %determines  the  number  of  initial 
%  codevectors 

c  =  initial  codebook; 


%Begin  the  algorithm 

D_avg_i  =  5 ; 

D_avg_i2  =  0; 

while  abs ( (D_avg_i-D_avg_i2 ) /D_avg_i)  >  epsilon  %Continue  while 

%  the  change  in  relative  error  is  still  greater  than  epsilon 

for  m  =  1:M  %Compare  each  image  to  each  codevector 

%  and  to  determine  which  codevector  the  image  belongs  to 

difference  =  zeros ( size ( snapshots , 1 ), size ( snapshots , 2 ), 1 ) ; 
min_dif  =  zeros (1,T); 
for~n  =  1:T 

dif ference ( : , : , n)  =  abs ( snapshots (:,:, m) -c (:,:, n) ) ; 
min_dif (n)  =  ... 

dot (reshape (difference ( :  ,  : ^  n)  ,  [ ] a  1 ) a  reshape (difference ( : ,  : , n) ,  [  ]  ,  1 )  )  ; 
end 

approx ( : , : , m)  =  c  \ ,  find (min_dif==min (min_dif ) ) ) ;  %build  an 
%  image  set  with  the  selected  codevectors  approximating  snapshots 
which_c (m)  =  f ind (min_dif==min (min_dif ) ) ; 

end 

for  n  =  1:T 

c ( : , : , n)  =  ... 

sum (snapshots  ( : ,  : ,  find (which_c==n) ) , 3) / size (find (which_c==n) , 2) ; 
end 

D_avg_i2  =  D_avg_i; 

D_avg_i  =  0 ; 

for  p  =  1:M  %Find  the  error  associated  with  the 

%  approximated  image  set 
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D_avg_i  =  D_avg_i+dot ( reshape (( snapshots (:,:, M)  ... 

-approx ( : ,  : , m) ) , 1 ,  [] )  r reshape ( ( snapshots ( : ,  : ^ M) -approx ( : ,  : , m) ) , 1 ,  [ ] ) ) 
end 

D_avg_i  =  1/M/ ( size  ( snapshots ,  1 )  "^size  ( snapshots ,  2 )) '^D_avg_i ; 

end 

basis  =  c; 
which  =  which  c; 
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Bestratio: 

function  [bestratio,  ratiostotry, sumresd] =bestratio (R, C_list,  . . . 
density, T, rratio, maxT, params , data, Beams ) 

%  function  to  determine  the  amount  by  which  to  scale  the  cfd  radius 
%  It  returns  the  ratio  which  had  the  smallest  residual  when  fit  to  the 
%  data. 

o 

o 

%  Author:  Ryan  Givens 

ratiostotry= ( 75 : 125 ) ./lOO; 
sumresd=ones (1, length ( ratiostotry)  )  ; 

for  i=l : length ( ratiostotry) 
i 

[A,  B,  X,  y]  =burnerbasisv2  (C_list ,  density,  T,  rratio,  params ,  R"^  ... 
ratiostotry ( i ) ,maxT) ; 
m=mean (B, 3 ) ; 

[basis , eig] =Create_basisv2 (B, m, 2 ) ; 

[ f it , resd, err ] =dataf it2 (m, basis , data ( : , 1 ) , x, y, Beams , data (:,2)); 
sumresd ( i ) =sum ( resd . ^2 ) ; 

end 

figure  (8009) 

plot (ratiostotry, sumresd) 
b=find ( sumresd==min ( sumresd) ) ; 
bestratio=ratiostotry (b) ; 
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Densmap: 

function  [nmap, nerrmap] =densmap (fit, fiterr, tmap, terrmap) 

%  function  to  create  a  density  map  from  the  temperature  map.  Errors 
%  are  propogated  to  the  nerrmap 

o 

o 

%  Author:  Ryan  Givens 

Sref=l . 168e-20; 
nu=7161 .4101; 
k=. 695036; 

E0=224 . 838; 
nmap=ones (size (fit) ) ; 
ner rmap=nmap ; 

for  i=l : size  ( fit , 1 ) 

for  j =1 : size  ( f it , 2 ) 

S=Sref ^partition (300, nu) /partition (tmap (i, j ) , nu) ^  ... 
exp (-EO/k/tmap (i, j ) ) /exp (-EO/k/300) ^ (1-exp (-  ... 
nu/k/tmap (i, j ) ) ) / (1-exp (-nu/k/300) ) ; 
nmap ( i , j ) =f it ( i ,  j ) /S ; 

nerrmap  ( i ,  j  )  =abs  (fiterr(i,j)  /S-fit(i,j)  (Sref"^  .  .  . 

partition  ( 300 ,  nu)  /exp  ( -EO/k/300 )  /  ( 1-exp  ( -nu/k/300 )  )  . 

exp (-EO/k/tmap (i, j ) ) /partition (tmap (i, j ) ,nu) ) ^- . . . 

2 ^partition (tmap (i, j ) , nu) ^exp (-EO/k/tmap (i, j ) ) ^ . . . 

EO/k/  (tmap  (i,  j  )  )  ^2/  (partition  (tmap  (i,j)  ^nu)  )  ^2'^Sref'^.  .  . 
partition (300, nu) / exp (-EO/k/300) / (1-exp (-  ... 
nu/k/300) ) ^terrmap (i, j ) ) ; 

end 

end 
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Tophatfixedcold: 

function 

[basis, resdmap, x, y] =tophat (centers, diameters , par ams , datal397 , . . . 

Beams , sigma) 

%function  to  create  an  optimized  two  zone  basis  image  with  the  cold 
%region  fit  only  to  beam  paths  not  passing  through  the  flame. 

o 

o 

%input : 

%centers :  array  of  center  locations  about  0  to  try 
%diameters :  array  of  diameters  to  try 
%params :  [Lx, Ly, dx, dy] 

%al397:  1397  integrated  absorptions 
%Beams :  Beam  positions 

%sigma:  errors  in  integrated  absorptions 

o 

o 

%Output : 

%resdmap:  2D  map  of  sum (abs (resdiduals) ) 

o 

o 

%Author:  Ryan  Givens 

resdmap=ones (length (centers) -1, length (diameters) ) ; 
normpathlength=sqrt ( (Beams (1,1,1) -Beams (1,2,1))^2  ... 

+  (Beams (2,1,1) -Beams (2,2,1) ) ^2 )  ; 
for  i=l : length (centers) -1 
i 

for  j =1 : length (diameters ) 

[basisa, x, y ] ^burner ( [ centers ( i ) , centers ( length ( centers ))],... 
diameters ( j ) , 2 , params ) ; 

[basisb, X, y] ^burner ( [centers (i) , centers (length (centers) )],... 

diameters ( j ) , 1 , params ) ; 
basisc=basisa-basisb; 

basisc ( : , : , 2) =basisa-2 ^basisc (:,:,!); 

%get  cold  basis  function 

basiscold=basisc ( :  ,  : , 2) ^ (datal397 (1) +datal397 (7) ) /2 /normpathlength; 

%get  projection  matrix  (integrated  area  of  cold  path  length  of 

each 

%path) 

[cond, projmat ] =Cond_Proj  (Beams , basiscold, x, y, 300 ) ; 

%subtract  out  cold  contribution  to  integrated  absorptions 
temp_abs=datal397-pro jmat ; 

%get  hot  basis  function 
basishot=basisc (:,:,!); 

%fit  to  hot  basis  function 
[bestf  it ,  resd,  f  item  ]  =  ... 

dataf it3 (basishot , temp_abs , x, y, Beams , sigma) ; 

%bestf its_arr ( : , : , i) =bestf it+basiscold; 


o 

o 

[basisa, x, y] ^burner ( [centers (i) , centers (length (centers) ) ]  ,  diameters ( j  )  , 
2 , params ) ; 

o 

o 

[basisb, x, y] ^burner ( [centers (i) , centers (length (centers) ) ] , diameters ( j ) , 
1 , params ) ; 

%  basis=basisa-basisb;  %top 
%  basis  (  : ,  : , 2) =basisa-2 ^basis (:,:,!);  %bottom 
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[bestfit^resd, fiterr] =dataf it3 (basis, datal397 , x, y, Beams, sigma) 
resdmap (i, j ) =sum (sum (abs (resd) ) ) ; 

end 


end 

x=diameters ; 

y=centers (1 : length (centers) -1)  ; 

[centerf it, diameterf it] =f ind ( resdmap==min (min (resdmap) ) ) ; 
bestcenter= [centers (centerfit) , centers (length (centers) ) ] 
bestdiameter=diameters (diameterf it ) 

[basisa, x, y ] ^burner (best center, best diameter, 2 , params ) ; 
[basisb, x, y ] ^burner (best center, best diameter, 1 , params ) ; 
basis=basisa-basisb; 
basis ( : , : , 2) =basisa-2 ^basis (:,:,!); 
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